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,  i'Measuremencs  of  cross-stream  and  streamwise  vorticity  have  been  made  over  important 
regions  of  an  experimental  simulation  of  the  bursting  process  of  turbulent  boundary  layers. 
Vorticity  measurements  have  also  been  made  of  the  starting  vortex  of  an  airfoil,  in  a  vortex 
TLog,  and  in  a  Stokes'  layer.  Using  a  new  technique  developed  under  this  contract,  these 
aeasurements  have  been  made  at  approximately  fifty  simultaneous  positions  over  the  flow 
structures  of  interest.  The  technique  can  accurately  measure  vorticity  and  strain  rate,  as 
instantaneous  Reynolds  stress  and  velocities  over  a  field.  It  requires  only  a  clock 
for  calibration.  -Furthetiaore,'  it  Sffectly  measures  these  quantities,  avoiding 
n  irect  interpretations.  It /employs  long  time  persistence  of  irradiation  of  a  photo- 

nark  fluid  particles.  The  technique  has  been  named  LIPA,  an,.acronym  for  Laser 
uced  Photochemical  Anemometry.,  Its  accuracy  has  been  measured  by  comparing  it  to  an 
cc  solution  of  the  Navier-Stokes  equations.  This  indicated  an  absolute  accuracy  of  cross 
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19.  ABSTRACT  (Continued) 


•rstream  vorticity  of  f^/sec.  Data  reduction  procedures  have  been  developed 
which  utilizeyhigh  resolution  digitization  and  image  processing.  Algorithms 
have  been  written  and  tested  which  allow  almost  complete  automation  of  the 
data  reduction  procedure. 
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INTRODUCTION 

Th«re  is  a  need  for  measurements  of  flow  quantities  such  as  vortlcity  and 
strain  rate  as  well  as  instantaneous  Reynolds  stress  over  and  area  in  turbulent 
and  unsteady  flows.  Examples  include  all  aspects  turbulence  structure  and 
entrainment  studies,  external  and  internal  aerodynamics,  studies  of  turbulence 
control,  of  mixing  and  combustion.  In  many  cases  there  is  a  need  to  obtain  more 
detailed  information  about  these  flows,  but  our  understanding  of  a  large  number 
of  simpler  flows  would  also  be  Increased.  Existing  probe  techniques  are  at  best 
able  to  provide  single  point  estimates  of  gradient  information.  Of  particular 
importance  is  the  need  to  develop  instinjusentation  to  extend  these  capabilities  to 
high  speed  flow  studies.  Another  need  is  for  instrumentation  techniques  that  can 
enable  quantitative  Lagrangian  information  to  be  obtained.  Evidence  of  these 
needs  is  the  very  high  level  of  activity  both  in  the  US  and  in  Europe  aimed  at 
developing  quantitative  field  measurement  techniques.  It  is  well  known  that 
measurements  of  vortlcity  using  a  hot-wire  probe  at  a  single  location  are  very 
difficult,  but  significantly  less  difficult  than  use  of  LDA  to  perform  the  same 
task.  Thus,  the  goal  of  obtaining  a  lar^e  number  of  simultaneous  vortlcity 
measurements  over  an  area  (that  ma..y  include  an  evolving  coherent  motion,  for 
example)  requires  a  different  approach. 

The  technological  transfer  between  physicists,  plasma  physicists, 
chemists,  computer  scientists  and  fluid  dynamicists,  has  been  the  basis  for  the 
application  of  new  techniques  to  measure  these  quantities  at  single  or  multiple 
points.  Several  new  techniques  for  turbulence  measurements  in  fluids  have  thus 
been  adapted  from  physics  and  chemistry.  They  are  based  on  atomic  or  molecular 
properties  of  the  medium  and  changes  in  these  properties.  The  techniques  have 
rapidly  become  useful  tools  because  of  advances  in  laser  and  computer  technology. 
Typical  examples  of  the  diagnostics  are:  Laser  Induced  Fluorescence  (LIF) ,  Forced 
Rayleigh  Scattering  (FRS).  Anti-Stokes  Raman  Diffusion  (ASRD) ,  Fourier 
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Densitometry,  and  Laser  Speckle.  These  advanced  techniques  can  lead  to 
Information  of  a  new  kind,  or  may  enable  measurements  to  be  performed  in 
currently  difficult  flows.  To  this  list  we  are  now  able  to  add  Laser  Induced 
Photochemical  Anemometry  (hereafter  call  LIPA) .  It  has  specific  advantages  over 
many  of  the  others,  which  are  discussed  below. 

Although  many  of  the  new  optical  techniques  are  essentially  single  point 
measurement  techniques,  LIPA  is  a  field  measurement  technique,  as  is  PV  (Particle 
Velocimetry,  which  includes  LPV,  Laser  Particle  Velocimetry) .  Of  those  above 
which  are  field  measurement  techniques,  all  require  optical  data  acquisition,  and 
thus  they  have  in  common  several  advantages  and  disadvantages .  Advantages  include 
either  no  intrusion  or  very  little  intrusion  on  the  flow  to  be  measured,  and  very 
fast  response.  Many  of  these  new  techniques  can  be  used  in  compressible, 
supersonic  or  hypersonic  flows,  where  the  conventional  techniques  (hot-wires, 
laser  doppler  anemometry)  are  difficult  to  use.  LIPA  can  also  be  used  in  high 
speed  flows,  although  these  applications  were  not  part  of  the  current  work. 

Common  disadvantages  include  the  need  for  optical  interrogation  of  the  results 
making  conversion  to  digital  form  expensive. 
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Host  of  these  new  field  techniques  suffer  from  limitations  of 
interpretation  (for  example  fourier  densitometry  provides  wave  vectors),  or  from 
limitations  to  low  speed  liquid  flows  (particle  velocimetry) ,  or  from  inherent 
inaccuracies  (for  example  the  spatial  ambiguities  of  speckle),  or  from  inaccurate 
data  reduction  potential  (holographic  particle  velocimetry) .  LIPA  can  surmount 
almost  all  of  these  difficulties.  The  outline  of  this  report  will  first  include  a 
discussion  of  the  technique  and  its  variations,  then  a  discussion  of  results  from 
the  present  investigation,  and  finally  an  overview  of  new  applications  and  the 
specific  advantages  and  disadvantages  of  using  LIPA  to  study  them. 
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RESULTS 

A)  The  LIFA  technique 

LIFA  depends  essentially  upon  Barking  specific  fluid  particles  so  that  they 
can  be  tracked  as  they  move  In  a  flow.  We  must  record  the  tagged  particles' 
motion.  We  have  restricted  our  data  acquisition  to  the  visible  range  of 
wavelengths,  but  the  technique  Is  not.  In  principle,  confined  to  this  range.  In 
the  visible  range  we  have  two  types  of  photochemical  tagging  available,  either 
color  change  or  emitted  light.  Color  change  upon  irradiation  is  called 
photochromlsm,  while  light  emission  is  broadly  called  florescence.  Both  types 
have  been  used  in  experiments  reported  herein,  and  thus  a  review  of  the  pros  and 
cons  of  each  are  now  presented. 

At  this  stage  in  the  development  of  LIFA|  a  limiting  factor  Is  the  ntunber  of 
chemicals  that  are  available  to  be  irradiated  with  lifetimes  long  enough  to  allow 
the  flow  to  move  appreciably.  This  limitation  has  made  the  lifetime  of  the 
chemical  govern  the  experiments  that  have  been  performed.  The  chemicals  have  also 
often  been  the  determining  factor  in  the  choice  of  fluid  that  has  been  used  for 
the  experiments.  In  general  photochromlc  chemicals  have  a  lifetime  of  Ims  to 
several  seconds.  This  lifetime  is  long  enough  so  that  they  are  suitable  for 
almost  all  experiments,  from  very  low  speed  to  very  high  speed.  However,  these 
chemicals  can  only  be  dissolved  in  organic  liquids.  Thus,  the  range  of  fluid 
dynamic  experiments  with  photochromlc  chemicals  is  limited  to  chose  in  liquids. 

On  the  ocher  hand,  photo luminescent  chemicals  tend  to  have  very  short  lifetimes. 
Typical  florescent  lifetimes  are  the  order  of  nanoseconds,  far  toe  short  for  the 
radiation  to  last  long  enough  for  all  but  hypersonic  flows  to  appreciably  move 
the  tagged  fluid  particles.  Thus,  special  photoluminescent  chemicals  are  needed 
which  have  longer  lifetimes.  These  tend  to  be  specially  constructed  molecules 
which  have  a  'shield'  of  ocher  atoms  surrounding  them  that  prevent  collisional 
quenching,  but  are  transparent  to  the  incident  and  emitted  radiation.  One  of 
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these  is  the  proprietary  product  called  Flowlite  I  from  Flowmod  Corp.  This 
chenical  can  be  dissolved  in  water  and  has  a  l.Sns  lifetime.  Ue  have  used  it  in 
preliminary  experiments  in  water  and  in  a  water  droplet  aerosoL  Another 
chemical  we  have  used  is  K2Pt2(POP)4-  This  was  synthesized  at  Michigan  State 
University  by  professor  D.  Nocera  It  has  a  lOps  lifetime,  and  has  been  used 
successfully  as  a  water  droplet  aerosol  in  N2.  Ve  envision  a  rapid  growth  of  the 
number  of  chemicals  synthesized  for  use  in  LIPA  as  the  value  of  the  technique 
becomes  apparent. 

For  the  experiments  described  herein,  tagging  has  been  performed  in  a 
plane.  This  is  not  an  inherent  constraint;  three-dimensional  tagging  and  stereo 
recording  is  certainly  possible,  and  is  seen  as  an  evolutionary  step  in  the  use 
of  LIPA.  Tagging  in  a  plane,  combined  with  stereo  viewing  will  enable  all  three 
components  of  velocity  in  the  neighborhood  of  a  chosen  plane  to  be  obtained. 
Furthermore,  creation  of  two  closely  spaced  grids,  combined  with  stereo  data 
acquisition  will  enable  all  three  components  of  the  vorticity  vector  in  the 
neighborhood  of  the  bisecting  plane  to  be  obtained. 

In  the  current  version  of  the  technique,  we  have  chosen  to  develop  passive 
beam  splitting  and  steering  devices,  which  divide  the  laser  beam  into  'n*  beams 
(we  have  used  n  -  7,  and  10)  and  cause  these  beams  to  cross  within  the  fluid.  It 
has  the  advantages  that  no  incident  light  is  lost,  and  that  the  line  widths  and 
spacings  can  be  specified  by  the  facet  design.  The  devices  are  similar  in 
appearance  to  an  oversized  diffraction  grating  (see  Appendix  B,  Fig.  1  and  2), 
but  it  gives  specular  reflection.  The  subsequent  excitation  of  the  chemical  along 
the  multiple  intersecting  laser  lines  creates  a  grid  of  'n'  lines  of  marked 
fluid.  At  each  of  the  intersections  of  the  laser  gridlines  a  'fluid  particle'  is 
unambiguously  tagged.  This  results  in  n^  simultaneous  measurements  over  the  area 
of  the  grid.  The  distortion  of  the  grid  in  a  known  time  provides  us  with  all  of 
the  kinematic  information  obtainable  in  a  fluid  In  the  plane  of  the  grid.  The 


resolution  of  this  information  is  equal  to  the  applied  grid  mesh.  For  flows  where 
differential  resolution  is  natural,  the  grid  can  be  set  to  match  the  needs  of  the 
flow. 

Details  of  the  use  of  the  technique  have  been  published  (see  Falco  and  Chu 
(1987);  Falco,  Chu,  Hetherington  and  Gendrich  (1988))  and  are  included  in 
appendix  B  and  C.  Further  description  is  given  in  Chu  (1988)  and  Falco  and  Chu 
(1988).  It  is  important  to  emphasize  that  the  technique  used  requires  only  a 
known  length  and  a  known  time  between  exciting  the  chemical  with  a  laser  and 
imaging  the  marked  fluid  particles.  Precision  placement  of  measurement  points  is 
another  advantage  of  the  technique.  Because  the  laser  excitation  is  in  a  plane, 
our  information  will  be  at  selected  points  in  (or  very  close  to)  selected  planes 
in  the  flow.  This  ability  to  choose  to  have  the  measurement  points  near  a  center 
line  or  at  y'*^  -  5,  for  example,  is  extremely  helpful. 

B)  Automated  data  reduction 

The  success  of  the  development  of  LIPA  into  a  cool  that  can  be  widely  used 
for  detailed  fluid  dynamic  measurements  will  most  probably  depend  upon  the  ease 
of  data  reduction.  Thus,  we  have  put  considerable  effort  into  automating  the  data 
reduction.  Key  to  our  success  has  been  the  use  of  the  lines  of  radiation  from  the 
chemicals  in  the  grids  to  provide  redundant  information  for  the  determination  of 
the  points  of  intersection  of  the  laser  lines.  This  has  made  imperfections  in  the 
film  or  optics,  film  grain,  dirt  acquired  in  developing,  etc.  unimportant.  The 
fact  that  two  independent  pieces  of  film  are  used  further  decouples  the  problems 
with  parasitic  images,  and  as  a  result  we  have  a  robust  data  reduction  technique. 

The  algorithms  used  are  included  in  the  Appendix  A.  They  are  written  in 
FORTRAN  and  have  been  run  on  a  Silicon  Graphics  IRIS  3120,  under  Unix  System  V. 
Briefly,  it  allows  the  velocity,  gradients,  vorticity  and  instantaneous  Reynolds 
stress  to  be  calculated  from  four  randomly  spaced  points.  The  field  may  have  'p' 


groups  of  points,  where  'p'  is  only  limited  by  available  computer  memory.  The 
algorithm  will  calculate  the  vorticity  in  contiguous  polygons,  if  points  are 
within  distances  that  are  less  than  a  specified  fluid  mechanically  determined 
limit,  such  as  a  Kolmogoroff  scale.  It  requires  that  each  of  four  points  be 
identified  in  two  photographs,  separated  by  an  arbitrary  time  Interval. 


C)  Determination  of  the  accuracy  of  the  technique 

Ve  have  calibrated  the  technique  by  making  measurements  in  a  flow  in  which 

we  have  an  exact  solution  (a  Stokes*  layer).  Results  are  presented  in  the  papers 

in  the  Appendix  B  and  C  (also  see  Falco  and  Chu  1988).  Essentially,  we  can  obtain 

cross-stream  vorticity  to  ±1  sec'^  accuracy.  This  compared  very  well  with 

estimates  made  using  classical  error  analysis  techniques.  Since  we  can  very 

accurately  measure  both  a  clock  and  a  ruler,  we  should  be  able  to  maintain  this 

accuracy  over  a  wide  range  of  flow  speeds. 


D)  Cross-stream  vorticity  measurements: 

Cross-stream  vorticity  measurements  have  been  made  in  the  Stokes'  layer  for 

calibration  as  mentioned  above,  in  a  vortex  ring,  and  in  the  starting  vortex  of 
an  airfoil.  The  Stokes'  layer  measurements  are  discussed  in  Appendix  B  and  C. 

The  measurements  in  a  vortex  ring  Falco  and  Chu  (1987),  show  for  the  first 
time  the  distribution  of  vorticity  in  a  single  ring  and  its  changes  as  we  move 
circumferentially.  Previous  measurements  with  LDA  have  required  averaging  over 
many  rings,  and  suggested  greater  symmetry.  These  changes  in  a  real  ring  that  was 
very  carefully  produced  and  very  stable,  indicate  the  insensitivity  of  the 
stability  of  the  ring  to  the  detailed  shape  of  the  vorticity  distribution  or  of 
its  circumferential  symmetry.  The  other  point  of  Interest  is  that  the  circulation 
for  our  low  Reynolds  number  ring  is  quite  close  to  that  found  for  a  Hill's 
spherical  vortex. 


The  measurements  of  the  starting  vortex  formed  when  an  airfoil  is 
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Inpulslvely  started  have  also  shown  significant  surprises.  This  study,  described 
In  Appendix  C,  Indicated  that  only  about  1/3  of  the  vortlclty  was  In  the  starting 
vortex,  which  means  that  the  remaining  vortlclty  Is  In  a  shear  layer  between  the 
wing  trailing  edge  and  the  vortex.  As  a  result,  the  Kutta  condition  can  not  be 
established.  This  has  Interesting  Implications  for  non-steady  airfoil  work  and 
for  calculations  that  depend  on  the  Kutta  condition  being  established.  Since  It 
has  not  been  possible  to  make  measurements  of  this  type  before,  unsteady 
phenomena  which  require  vortlclty  measurements  over  an  area  have  recieved  no 
experimental  attention. 

E)  Streamwlse  vortlclty  measurements  of  a  simulation  of  the  production  process 
In  turbulent  boundary  layers 

Streamwlse  vortlclty  measurements  have  been  made  In  streamwlse  vortices 
created  above  a  wall,  similar  to  those  found  In  turbulent  boundary  layers.  These 
results  have  been  reported  by  Chu  (1988)  and  by  Falco  and  Chu  (1988).  Briefly,  we 
have  observed  this  sequence  in  the  turbulent  boundary  layer,  and  have 
subsequently  been  able  to  model  it,  in  all  its  detail,  by  passing  a  vortex  ring 
over  a  laminar  shear  layer  (either  a  Stoke's  or  Blaslus  layer).  It  is  true  that 
the  complete  event  Is  rarely  seen  in  Che  turbulent  boundary  layer.  However,  a 
significant  part  of  It  was  apparent  during  each  of  the  bursts  we  have  observed 
(and  It  must  be  noted  that  our  techniques  don't  result  in  observation  of  the 
entire  picture  at  all  times).  We  are  presently  doing  this  in  the  turbulent 
boundary  layer  using  either  multiple  dye  silts,  or  the  LIPA  technique  (it  will  be 
reported  on  later) ,  which  Is  allowing  whole  field  real  time  measurements  of 
vortlclty  and  Instantaneous  Reynolds  stress  in  the  marked  fluid. 

In  experiments  designed  to  examine  only  the  wall  region,  we  observed  that 
the  production  sequence  starts  with  the  formation  of  a  pair  of  long  screaks. 

Using  two  dye  slits  in  series,  we  observed  that  subsequently  a  pocket  forms 


between  the  streak  pair;  often  after  streaks  of  several  hundred  x*  have  formed. 

It  appears  to  us  that  the  dye  marks  a  strong  vortex  which  forms  in  the  pocket. 

Use  of  the  LIPA  technique,  in  simulations,  has  clearly  showed  that  the  pocket 
vortex  is  present  and  that  it  is  the  most  intense  vortex  found  in  the  wall 
region. 

The  other  aspect  that  has  been  quantified,  is  that  the  long  streaks  are 
observed  to  undergo  an  instability  that  results  in  the  formation  of  one  or  more 
hairpin  vortices  over  each  streak.  Using  LIPA,  we  have  measured  the  streamwlse 
vorticity  in  these  vortices  and  found  it  to  be  weak  with  respect  to  the 
streamwlse  vorticity  in  the  pocket  vortex.  Left  to  evolve  on  their  own,  these 
hairpin  vortices  would  not  contribute  much  to  the  production  of  turbulence. 

However,  this  is  not  the  case.  The  experiment  shows  that  as  the  pocket 
vortex  is  stretched,  these  secondary  hairpins  are  convected  past  it.  Because  the 
pocket  vortex,  which  now  has  a  primarily  streamwlse  component,  is  much  stronger, 
it  induces  the  hairpin  vortices  around  it.  Our  two-dimensional  instantaneous 
vorticity  maps  of  these  features  and  their  interactions  gave  us  a  detailed 
quantitative  account  of  the  ensuing  interaction.  We  found  that  the  pocket  vortex 
was  more  than  twice  as  strong  as  the  hairpins,  and  that  it  lifted  them  completely 
around  it  leading  to  the  breakup.  Furthermore,  practically  all  of  the  Reynolds 
stress  associated  with  the  screak. . .secondary  hairpin. . .pocket  vortex  interaction 
is  due  to  the  pocket  vortex  and  the  motion  it  induces. 

In  addition,  our  observations  have  indicated  that  the  streaks  do  not  form 
as  a  result  of  a  'pumping'  action  by  long  streamwlse  vortices,  and  that  evidence 
of  vortices  associated  with  the  streaks  is  only  found  when  the  secondary  hairpins 
form.  LIPA  unequivocally  shows  no  concentration  of  streamwlse  vorticity  during 
the  time  when  screaks  first  appear. 


DISCUSSION 

There  are  two  additional  aspects  that  use  of  the  LIPA  technique  can 
uniquely  provide  for  researchers  in  turbulent  and  unsteady  flows.  They  are  1) 
quantitative  Lagrangian  Infomatlon  and  2)  quantification  of  the  kinematics  of 
high  speed  flows. 

Since  high  speed  movies  can  be  taken  (synchronized  in  combination  with  a 
high  repetition  rate  Excimer  laser)  we  can  gain  an  understanding  of  the  evolution 
of  the  dynamics  of  unsteady  fluctuating  flows,  i.e.  obtain  quantitative 
Lagrangian  information.  This  is  of  great  importance  in  the  study  o^  control  of 
turbulence,  for  the  details  of  how  the  control  interacts  with  the  coherent 
structures  is  the  essential  design  information.  Experiments  of  this  type  are 
currently  underway  at  the  Turbulence  Structure  Laboratory  to  help  understand 
turbulent  boundary  layer  inner-outer  interactions,  and  their  role  in  the 
production  process. 

Another  range  of  insights  will  come  from  far  more  detailed  measurements  of 
the  dynamics  of  high  speed  flows.  The  capability  of  LIPA  to  obtain  whole  field 
measurements  in  gas  flows  (which  are  seeded  with  submicron  water  droplets 
containing  one  of  the  chemicals),  makes  possible  measurements  in  high  speed 
flows.  Most  uv  lasers  pulse  in  nanoseconds,  and  the  reactions  occur  in  even 
shorter  time  scales.  Since  this  is  a  non- intrusive  measurement  technique,  and 
since  it  requires  only  a  clock  and  a  ruler  for  calibration,  it  potentially  can 
provide  high  accuracy  results  for  flows  at  any  speeds  currently  obtainable  in 
wind  tunnels.  Recording  the  signals  with  image  converter  cameras  appears  entirely 
possible.  Thus,  potentially,  transonic  and  supersonic  flow  studies  could  be 
performed.  An  interesting  sidelight  is  that  LIPA  may  prove  to  be  a  practical 
Cool  for  aircraft  developrent.  Since  many  high  Reynolds  number  experiments  are 
now  being  performed  in  cryogenic  tunnels,  and  since  Nj  essentially  prevents 
quenching,  construction  of  an  optical  setup  that  would  provide,  say  100  to  1000 
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data  points  over  a  body,  might  prove  very  cost  effect  when  compared  to  the  tunnel 
time  required  for  a  similar  number  of  LDA  measurements. 

Finally,  LIPA  using  a  gaseous  phase  chemical  in  a  gas  flow  (this  is 
conceptxially  possible,  but  has  not  been  attempted),  hes  the  potential  to  give 
whole  field  measurements  in  high  temperature,  combusting  or  hypersonic  flows. 

Suomiarizing,  the  advantages  of  the  LIPA  technique  are: 

•  Instantaneous  measurements  of  velocity  gradients,  vorticity, 

strain  rate.  Instantaneous  Reynolds  stresses  and 
velocities  can  be  made  over  an  area 

•  Can  be  setup  to  give  highly  accurate  space  and  time  resolved 

measurements  over  a  large  range  of  velocities 

•  Non- intrusive 

•  Measurements  can  be  made  at  chosen  locations 

•  Simple  calibration  depends  only  on  a  clock  and  a  ruler 

•  Simple  data  acquisition 

•  Not  limited  by  Interferometric  quality  access 

•  Data  reduction  can  be  automated 

•  Can  be  used  in  gases  (currently  in  aerosol  form)  or  liquids 

•  Can  be  used  in  non-Newtonian  fluids  and  in  two -phase  flows 

•  Three-dimensional  measurements  are  possible 

•  Standard  wavelengths  of  uv  lasers  can  be  used 

•  Can  be  combined  with  flow  visualization 

•  Mixing  can  be  studied  by  using  more  than  one  chemical 

•  Potentially  can  be  combined  with  LIF  to  give  concentration 

information  along  with  the  kinematic  data. 

The  chief  disadvantage  of  LIPA  is: 

•  Specific  chemicals  must  be  found  for  flow  velocities  and 

fluids  of  interest. 
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CONCLUSIONS  AND  FUTURE  WORK 

UFA  has  been  developed  and  used  to  measure  derivative  information  at 

multiple  points  in  complex  flows.  An  algorithm  has  been  written  that  essentially 
automates  data  reduction.  The  study  of  a  model  of  the  production  process  has 
shown  that  the  pocket  vortex  is  the  strongest,  easily  overwhelming  the  hairpin 
vortices  that  develop  over  the  streaks.  Our  study  of  the  starting  vortex 
Indicates  that  significant  vorticity  exists  in  a  shear  layer  between  it  and  the 
airfoil  at  low  Reynolds  numbers. 

Quantification  of  this  kind  has  here-to-fore  been  impossible.  LIPA  enables 
quantitative  measurements  of  velocity,  vorticity,  strain  rate,  and  instantaneous 
Reynolds  stress  over  an  area  in  turbulent  flows.  Only  two  successive  pictures  are 
needed  to  calculate  the  distortion  of  the  marked  fluid.  Calibration  requires 
only  a  clock  and  a  ruler.  It  works  in  reversed  flows.  The  technique  has  been  used 
in  kerosene,  water,  nitrogen  and  air  (with  seeding).  As  such,  it  should  be  able 
to  measure  flows  into  the  supersonic  range.  It  can  also  be  extended  so  as  to  make 
three-dimensional  measurements  of  all  kinematic  quantities  over  an  area.  The 
potential  for  the  use  of  LIPA  in  many  other  flows  and  particularly  at  high 
Reynolds  numbers  has  been  discussed. 

Future  work  will  include: 

a)  Measurement  of  coherent  motions;  their  formation, 

evolutions  and  interactions  in  turbulent  water  flows 

b)  Development  of  technique  in  gases. 

c)  Use  of  the  technique  in  water  to  develop  its 

three-dimensional  capability  for  both  velocity  and 
vorticity  fields. 

d)  Development  of  simultaneous  LIPA  and  LIF  measurements. 

e)  Use  of  the  technique  in  high  speed  flows. 
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This  appendix  contains  the  comiaann  files  which  run  the  image 
processing  equipment,  converting  che  visual  information  into  uigitai 
lont.  The  program  wnich  interrogates  the  visual  information  and 
calculates  cne  fluid  kinematic  quantities  is  also  inciuued.  The 
programs  which  are  run  on  che  image  processor  to  perform  standard 
operations  like  averaging  and  edge  ennancement  are  too  lengthy  to 
include  acre. 

The  data  flow  la  as  follows.  Either  film  or  vioeo  tape  is  usee 
to  record  an  experiment.  The  image  processor  is  requireu  to  digitize 
ail  of  tnese  frames  ana  Chen  iijdiviau*»ily  reouce  u-ach  image  to  a  set 
u!  intersection  points.  A  vaodel  coumanc  file  can  be  applieu  to  eacn 
-;,.age  to  perform  this  >-ata  reuucticn,  out  since  our  image  processor  noes 
not  allow  one  coamano  file  to  call  anotr.ar,  „  large  conmanc  file  u.ust  oe 
ogneracea  to  aucocaicicaily  process  nainy  iraiaes  without  personal  i.TCer- 
veution. 

The  rirst  command  file  (TapeCrunch)  is  a  UNIX  shell  script 
which  generates  tnis  large  image  processing  command  file.  In  this 
exai.-pie  it  IS  setup  to  acquire  up  Co  dO  sequential  frames  iroii  a  vioeo 
cape,  ..tore  cneu.  and  then  process  them  us in^  ..  -...odei  command  cile. 

The  secy.no  ccmmar.c  file  •TapeCruncn.cccr.;  is  suitaole  fur  use  on  che 
ilefaaVision  1024X1.  image  processor.  AithoUj,a  up  to  4C  frames  icay  be 
prucessea  at  once,  \;e  show  how  only  five  are  done  in  the  interests  of 
orevicy.  The  third  cemmana  file  lauto.cxml  contains  the  commancs  to 
reouce  sne  uij,itiaea  images  to  two  secs  or  lutersecciuo  lines. 

The  program  VORTICITY  and  ail  ot  its  subroutines  comprise  che 
remainaer  ..i,  cnis  append;..  It  must  ce  run  after  the  intersections 
nave  oeen  uotaineo.  It  cakes  the  output  ot  the  image  processing 
programs  ar.o  converts  c.cat  into  velocity,  velocity  derivative,  anu 
.orticicy  information. 
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u!  /bin/sh 

G 

Li  Description: 

0  This  Unix  shell  script  generates  a  large  MegaVision  command  file 

G  to  crunch  many  frames  from  a  video  tape.  It  is  necessary  to 

□  generate  one  large  MegaVision  command  file  since  the  image 

Q  processing  software  doesn't  permit  nesting  from  one  command 

Li  to  another  then  back  again. 

U 

G  Use  of  Arguments  $0.  $1.  $2t  ...: 

G  ~  The  command  file  (3F1LE}  created  will  iiave  the  name  $0.cxm. 

G  NOTE:  Any  previous  contents  of  $FILE  are  wiped  out  by  this 

□  shell  script. 

Li  -  The  root  filename  used  to  store  the  video  data  frames  is  cake  from 
G  SI.  If  $1  IS  not  set  (or  null),  the  string  'VidFram'  is  used. 

^  -  The  root  filename  used  to  score  tne  incerseccion  daca  will  be  $2  or 

i_i  'Intersect'. 

-  The  nooel  conmana  file  co  cruncn  the  data  will  be  $3  or  'auto.omi'. 
lj  History: 

u  <cpg>  6  oct  88  —  1st  version.  Always  did  it  by  nano  before... 
u 

FILE=$0,c:cm 
FRAMES =5 

arg=§l  ,  VROOT=‘>pwd®/${arg:=VidFram} 
arg=32  ,  IROOT=  »pwd<>/$ {arg :  =Intersect } 
arg=$3  ,  COMND=S{arg:=auto.cxm} 

ecno  ! 

!  Viaeo  tape  crunch  commanc  file. 

!  Read  in  SFRAMES  frames,  store  them,  and  then  recall  and  crunch 
!  eaca  lauividually  usiUj,  sCOtll.D  as  the  comnanu  file  template. 

I 

!  This  file  was  created  on  •=date®. 

I 

!  First  read  the  video  tape 

ERALL 

IHIT 

TVPACK  t  525,$FRAMES,0 

; 

!  Tiien  score  the  oata"  >  3FILE 

^  He  need  to  Ioo<c  at  Al ,  A2,  A3,  A4,  51,  ....  2  |  I 

until  we've  stereo  as  uany  frames  as  were  - + - 

L,  scannea  in.  The  frames  are  stored  like  this:  4  |  3 

ui  (See  tne  MegaVision  TVPACK  help.)  This 

G  uetermines  the  values  for  Xoff  ana  Yoff  and  the  order  in  whicn 
LJ  tney  change... 

F?d'(  _i<UK=0 

tor  ;;EM0P.Y  m  Al  A2  A3  A4  Ll  22  B3  34  C  D  E  F 
uo 

for  Yoff  in  0  512 
00 
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for  Xoff  in  512  0 
do 

FRM_NUM=oexpr  $FR1I_HUM  +  1“ 
if  test  $FR11_NUM  -gt  $FRAMES 
then 

break  3 
f  i 

FMAME= $VR00T$FRM_NUW 

echo  “USTORE  §MEM0RY  $FNAME  1 .312,512.$Xoff .$Yoff .1.1”  »  SFILE 
done 
done 
none 

echo  ”! 

!  DATA  STORED. 

!'■  »  $FILE 

FRlI_llUl-i=i 

while  test  SFRl-UIUW  -le  SFRAl-IES 

bO 

VKAME=  $VR00T$FR11_1}UM 
INAME=  iilROOT .  SFaK_i«JM 
echo  ”! 

!  Now  processing  frame  number  $FR11_NUM 

RECALL  SVNAIIE  B  0.0. 2. 2 
AVERAGE  B  A4  3. 5. 0.1 
scat  SCOllKD® 

I 

!  Calculate  and  store  the  intersections 
^OTS  B./id  C.Sli^Al'IE 
!"  »  SFILL 

FRli_NUM=<'expr  SFRIU'IUM  -  1» 
done 
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!  Video  Cape  crunch  cocrsand  file. 

!  Read  in  3  frames*  store  them,  and  then  recall  and  crunch 
!  each  individually  usin^'  model. cxm  as  the  command  file  template. 

t 

!  This  file  was  created  on  Thu  Oct  6  11:01:36  EOT  1988. 

I 

!  First  read  Che  video  tape 

ERALL 
INIT 

TVPACK  t  525,5,0 

I 

!  Then  store  the  data 

NSTORE  Al  /usr/thor/gendrich/vis j on/Videoi 
USTORE  Al  /usr/thor/gendrich/vision/Video2 
dSTORE  Al  /usr/thor/gendrich/vision/Video3 
NSTORE  Al  /usr/thor/gendrich/vision/Video4 
liSTORE  A2  /usr/ thor/'genarich/vision/Viaeo5 
! 

!  DATA  STORED. 


1.512.512.512,0,1,1 

1.512.512.0.0.1,1 

1.512.512,512.512,1,1 

1,512.512.0.512.1.1 

1.512,512.512.0,1.1 


!  Wow  processing  frame  number  1 

I 

RECALL  /usr/chor/genarich/vision/Viueol  B  0,0, 2, 2 
AVERAGE  3  A4  3, 5. 0,1 


.  model  command  file  is  inserted  here 
.  processing  aata  from  memory  A4  and 
.  uepositing  the  intersections  in  Bl . 


!  Calculate  ana  score  the  intersections 
DOTS  5,A4  C,/usr/tnor/gendrich/vision/Intersecc . 1 
; 

I 

!  liow  processing  frame  number  2 

j 

RECALL  /usr/thor/gencricn/vision/Video2  B  0,0, 2, 2 
AVERAGE  B  A4  3, 5. 0.1 

•  ^rodei  commanu  file  is  inserted  here 
.  processing  data  from  memory  A4  and 
.  depositing  the  intersections  in  El . 

I 

!  Calculate  and  store  the  intersections 

DOTS  D,A4  C./usr/thor/genorich/vision/Intersect.E 

I 

; 

!  Wow  processing  frame  number  3 

I 

RECALL  /usr/ thor/gendrich/vision/Viaeo3  3  0,0, 2, 2 
AVERAGE  L  A4  3, 5, 0,1 
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.  model  command  file  is  inserted  here 
.  processing  data  from  memory  A4  and 
.  depositing  the  intersections  in  El. 

• 

I 

1  !  Calculate  and  store  the  intersections 

DOTS  B,A4  C,/usr/thor/gendrich/vision/Intersecc.3 

I 

! 

!  IJow  processing  frame  number  4 

I 

b  RECALL  /usr/thor/gendrich/vision/Viaeo4  D  0.0. 2<2 

r  AVERAGE  B  A4  3. 5. 0.1 


► 


.  model  command  file  is  inserteo  iiere 
.  processing  data  from  memory  A4  ann 
.  depositing  the  intersections  in  El. 

« 

J 

!  Calculate  and  store  the  intersections 

DOTS  B,A4  C,/usr/thor/gendrich/vision/Intersect.4 

I 

! 

!  Now  processing  frame  number  5 

I 

RECALL  /usr/chor/gendrich/vision/VideoS  E  0.0. 2. 2 
AVER/.GE  B  A4  3. 5, 0,1 


;:oaei  commanc  file  is  laserteu  iiere 
processing  aata  from  memory  a4  ana 
aepositing  tne  intersections  in  Bl . 


Calculate  and  store  the  intersections 
OTS  E,A4  C,/usr/chor/gendrich/ vision/Intersect . 5 
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luage  processing  conmana  file 

Works  with  all; 

Ektachrome  wing  images 
ASA  1000  wing  images 

NOTE;  The  angles  through  which  the  picture  is  rotatea  are  dependent  on 
the  setup  of  the  beam  dividers.  These  values  will  change  from  one 
setup  to  the  next,  requiring  modification  of  the  'ROTATE'  angle  at 
the  start  and  end  of  processing  each  half  frame. 

History; 

<cpg>  1  Jan  88  —  comments 

<cpg>  21  Jan  88  —  Instead  of  subtracting  the  oackground 
level,  differ  ana  HEQ  are  used. 

PRELIMINARIES 


!  first,  subtract  the  backgrouno  level  to  get  a  more  even 

!  grey  level  across  the  entire  screen. 

average  a4  b4  1,37,3,0,1 

subtract  34,04  c 

threshold  c  a3  0,143,2,2,2,2 

I 

!  To  eliminate  the  noise  at  the  edges,  crop  A3, 
erase  b3 

block  a3  b3  64,64.959,959,64,64 

move  b3  a3 

I 

! 

!  START  PROCESSING  THE  FIRST  ilALF  FRAME; 

!  Rotate  the  image  so  the  first  set  of  lines  is  vertical 
rotate  j3  b3  -25 

I 

!  Eapnasize  the  vertical  lines  --  their  average  within  a  vertical 
!  kernel  is  4,reater  than  the  surroundings  (and  the  other  lines), 
average  b3  a2  1,1,130,0,1 

tnresnold  a2  b2  0,168,2,2,2,2 

not  b2  al 

t 

!  Fill  in  gaps  from  thresholding. 

_range  al  bl  1,1,256,100,255,255 

I 

!  Eliminate  tne  thin  noise  spikes, 
minimum  bl  a3  3,5 

I 

!  Fatten  up  the  lines  a  bit  --  it  helps  for  filling  in  the  gaps, 
rollor  a3  53  3 

[ 

!  Now  fill  in  tne  gaps. 

_range  o3  a2  1,1,256,110,255.255 

rollor  a2  52  5 

I 

J  Elirrinate  any  stray  hign  order  bits... 
tnresnold  52  ^  . 
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!  Shrink  lines  to  one  pixel  widet  being  careful  to  eliminate  any 

!  open  circles  and  line  stubs.  (ERODEC  erodes  everything  but  Contiguous 

!  objects.) 

move  al  bl 

erodec  bl  a3  20 

!  (This  fills  in  holes  for  any  points  that  might  not  want  to  go  away...} 
rollor  a3  b3  3 
!  Now  finish  up  the  job... 
erodec  b3  a3  36 

I 

!  Put  these  lines  into  their  original  frame  of  reference, 
rotate  a3  b3  25 

j 

!  t'.OTE  —  the  origin  wasn’t  coved  rn  our  croppiu^  above  (although  it 
!  could  have  been).  If  the  origin  *is*  moved,  use  CLOCK  here 

!  to  put  it  back  where  it  belongs. 

!  1  2 


Next,  fatten  up  tne  lines  so  that  intersections  1  2 

won't  be  missed  by  one  line  zigging  while  1  2 

the  other  is  tagging.  See  illustration  ==>  12 

21 

The  intersection  between  lines  "1"'  ana  2  1 

"2”  IS  a  <NULL>  instead  of  a  point.  2  1 


rollor  b3  a2  3 

END  OF  PROCESSING  FIRST  liALF  FRAiME. 


!  ...  now 

CO 

it 

all 

over  with  the  lines 

in  the 

other  direction 

thresnolu 

a2 

b2 

0.150.2,2.2,2 

not 

b2 

al 

suotract 

bd. 

a4 

C 

uiove 

1 

al 

b4 

store  the  tirst 

licilf  ir 

ame  in  24  wneu  taac 

[ 

raemory  is  aval  It 

ible 

cr.resnola 

c 

a3 

0,123.2,2,2,2 

erase 

b3 

olock 

a3 

b3 

64,64,959.959,64.64 

move 

b3 

a3 

I 

!  STAF.T  OF  PROCESSING  THE  2no  ilALF  FRAIIE 


rotate 

a3 

b3 

27 

average 

13 

a2 

1,1.100.0,1 

tnreshold 

a2 

b2 

0,166,2,2,2,2 

not 

b2 

al 

-range 

al 

bl 

1,1.256,110,255.255 

minimum 

bl 

a3 

3,5 

rollor 

a3 

b3 

3 

_range 

b3 

a2 

1.1,256,110,255,255 

rollor 

a2 

b2 

7 

tnresnold 

b2 

al 

0,150,2,2,2.2 

not 

al 

bl 

erodec 

d1 

a3 

20 

rollor 

a3 

b3 

3 

erodec 

o3 

a3 

36 

rotate 

a3 

b3 

-27 

rollor 

o3 

a2 

3 

Oct 
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!  END  OF  PROCESSING  THE  2nd  HALF  FRAME 

j 

!  Here  are  the  other  lines, 
tnreshold  a2  b2  0»1S0.2«2.1 ,2 
not  b2  al 

t 

!  A  logical  AND  is  the  intersection  ot  the  two  sets  of  lines 
!  (literally).  :-) 
ana  al«b4  c 

I 

!  The  system  uses  C  as  a  working  memory.  Put  our  intersections 
!  somewnere  safer.  :-| 
move  c  b  1 

I 

!  Let’s  overlay  the  points  on  the  original  image. 

!  VJhite  oojects  in  the  graphics  memory  \.F)  are  unat 
!  we  need  to  ao  that, 
i.iove  c  f 

I 

Look  at  tne  original  mage... 
view  a4 

I 

!  ...ana  turn  on  the  grapnics  display. 

oOn 
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I 


I 


> 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


program  VORTICITY 
Description: 

Ibis  program  reads  two  data  files  produced  by  the 
‘object*’  function  of  the  MegaVisiun  1024X11.  The  labelled  ana 
filtered  objects  are  tne  points  within  the  flow  which  we  are 
tracking. 

This  program  performs  a  number  of  tasks: 

1)  reads  the  data 

2)  sorts  it  by  increasing  X  and  Y  values 

31  matches  points  in  frame  1  to  points  in  frame  2 
4;  interpolates  iu  space  and  time  to  construct  a  velocity  field 
throughout  the  frame 

3)  connects  the  dots  to  form  polygons. 

6;  uses  the  velocity  inforoation  to  oetermine  derivatives, 
stress  and  strain  values,  vorticity.  etc. 


c  Autnur:  Chuck  Gendricn  i.Give  (r.e  unix  or  .../ 


c  History: 

c  August,  1987  --  1st  version  --  written  on  freyja.  an  IRIS  312C 

c  running  unix  System  V. 

c  Decemoer,  b7  --  1.1  <cpg>  Acdea  i,rapnics  display  or  the 

c  calculations  ^velocities,  vorticity,  etc.; 

c  Added  TURB3D-f ormatted  output, 

c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

CD  aefine  DEBUG 

include  "vorticity. n* 
c 


.ogicai  AoUin  .  I'uDiu: 
integer  .'umPts,  .iumPoly 
external  Again,  liumPts,  CumPcIy 

c  Again  is  the  error  hanuler  --  prints  oiagnostics  ana  asxs 

0  if  the  user  wants  to  quit  or  re-run  the  preceaing  section. 


c  UnDim  comes  back  TRUE  wnen  output  will  be  non-dim'd 

c 

openv  unit=l,  f ile- '/dev/tty • ) 

c  unitl  —  tty  input 

open;  unit  =  2,  file= '/dev/tty • ; 

c  unit2  --  tty  output 

c 

ccccccccclccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc 
c  step  1  —  get  Che  data 


u-  i:iciude  “orapnics.h” 


iC  continue 

call  getfrmv  iramel,  frames.  Uci’im.  lerr; 
if;  lerr.ne.O;  then 
c  error  nanale 

if;  Again;  lerr;)  goto  10 


enaif 
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c 

c 

c 

c 


J  itdet 


1000 


Li  enaif 


c 

c 

c 


20 


u  ifdet 

1010 
Ij  endif 
c 
c 
c 

L-  ifcer 

1020 

+  ' 
T 

enoif 

c 

c 

c 

30 


u  ifaet 
1030 


1040 
u  enaif 
c 


step  2  —  sort  the  points,  this  is  necessary  when  using  more 
sophisticated  matching  algorithms. 


call  sortfrmv  frarael) 
call  sort£rm(  £rame2) 

DEBUG 

u  =  NumPtsC  £ramel.  2) 
writer  *1  10001  'frame  l',n> 
n  =  NumPtsC  frame2.  2) 
writeC  *>  1000;  'frame  2',n. 
format (/'  X-Y  data  for  '.a7, 
',i2.':  ('.£8.3,', ',£8.3, 


\i,framel(i,X),framel(i,Y) ,i=l,n; 

( i .£ rame2( i .X) .frame2C i .Y } , i=l ,n / 
; ' ,i2 , '  points ; : ’/ 

)')) 


step  3  —  define  the  correspondence  between  the  two  frames  (match 
cont inue 

call  color!  GREEN! 

call  MatcnesC  framel,  frame2,  .natch,  terr; 
ifC  lerr.ne.Ol  tnen 

i£(  Again!  terr!)  ^oto  20 

end  if 
DEBUG 

writeC  1010}  match 

formate '  Here  is  the  match  array ;'/( i3 ) ) 


step  4  —  interpolate  the  velocity  vectors 

call  Velocities^  framel,  irameZ,  match,  vel} 

DEBUG 

11  -  .'iuciPcsv  vet,  4; 

writei  *,  1020;  . i , vet v  i ,X} . vel! i.Y )  ,vei!  i,  Vx; , vel!  i ,  Vy ; ,  i  =  l  .i. , 
formate/'  Here  are  the  velocity  components  and  their', 
locations;'/'  X  Y  Vx  Vy'/ 

(lx,i2,'.  ' ,4eix,f7 .2} } } 


step  5  —  define  the  polygons 
cont inue 

call  Polygons!  vel,  poly,  ierrj 
if!  lerr.ne.O;  then 

ifC  Again!  lerrj)  goto  30 

»'ndif 

DL?UG 

write!  *,  1030) 

format!/'  Here  are  the  polygons  wnich  have  been  defined;'// 
ULH  URH  LRH  LLH'} 
n  =  NumPoly!  poiy) 

writet  *,  1040,  .polyt i,l } ,polyi i,2} ,poiy! i,3; ,poly! i,4) , 
1-1 .n ; 

forma t  !4t  2:t,  i2 ,lx  / ; 
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L  ioclude  "Vusr/iiiclude/fgl.h” 

Li  include  "/usr/include/f device. h" 
c 

c  data  declarations  for  the  vorticity  program 

c 

integer  MaxPts,  flaxPoly 
parameter(  IlaxPts  =  100) 

c  the  maximum  number  of  points  per  frame 

c 

parameter C  llaxPoly  =  MaxPts; 
c  Che  maximum  number  of  polygons 

c 

real  LIGRADi  SIlALLIuiD 
parameter (  BIGRAD  =  3.0/ 
parameter C  SMALLRAD  =  1.5) 

o  large  and  small  radii  for  outer-  ana  inner-circle  crawii.g 

c 

real  4,4/,  sca..e 

a  cranscornation  tmcrix  iCr  sizing  things  correctly 

a  and  the  value  to  put  on  tne  uiagcnai 

c 

integer  X,  V 
parameter (  X  -  1/ 
paramecerv  7-2) 

c  pointers  to  XY  coords 

c 

integer  Vx,  Vy 
parameter!  Vx  =  3) 
parameter!  Vy  =  4) 

c  pointers  to  velocity  components 


c 

c 


c 

c 

c 

c 


c 

c 

c 


integer  BAD 
parameter!  BAD  -  0; 

wnen  .oox'  returns  a  BAD  point,  it  cian’t  fina  u 
corner  in  the  indicoteu  uirection  at  the  correct 
uiscance . 

integer  ddx,  ddy 
t/arameter!  ddx  =  X) 
parameter!  udy  =  Y; 

pointers  so  that  subprograms  agree  wnecher  the 
derivative  wrt  X  ur  wtc  Y  should  oe  cakei. 


real 

iraue! 

.iaxPcs , 

1) 

real 

irauel ! 

liaxPts . 

2) 

real 

.rameZ! 

MaxPts , 

2/ 

XY  coords  of  each  point  in  frames  1  ana  2 

integer  rratchC  haxPts; 

framelvi/  <-->  frame2(  match! i/l 

integer  poly!  XaxPoly,  4; 

making  four-siaed  polynomials 

poly!  1,  i,  is  tne  upper-Lli  corner  ci  tr.e  polyoou 
poly!  1,  2)  is  the  upper-kH  corner 
poly!  1,  3/  lowerLil 
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c 

c 

c 

c 


c 

c 

c 

c 

c 


c 

c 


0 

c 

c 

c 


poiyv  i.  4)  lowerRH 

real  vei(  MaxPts,  4) 

velocicy  tield  -  Vx  a  Vy  C  X,Y 

integer  I-latchDim 
parameter (  HatchDim  =  10; 
integer  matchs(  HatchDim; 

watch  the  spelling!!! 

Hopefully,  no  uore  than  10  points  will  match... 

characterise  title 

for  titling  output 

real  LP,  DOWN,  RIGHT,  LEFT 
parameteri  UP  =  90.0; 

parameter!  DOWN  -  270.0; 
parameter!  RIGHT  -  0.0. 

parameter!  LEFT  -  lOO.Uy 

integer  Xsize,  Ysixe 

size  of  tae  window  port  cpeneci  up. 

logical  GoOn 
external  GoOn 

a  logical  ^unction  which  types  ! <RET>  to  continue;  then 
returns  .TRUE. 

integer  ierr 

returns  0  v;uen  there  is  no  error,  otherwise... 

integer  i.  j,  k,  n,  nl ,  a2,  p 
counters 


tcccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


Apr  7  10; 13  1988  graphics. n  Page  i 


c 

c 


c 

c 

c 

c 

c 

c 

c 


77777 


c* 

c* 

c*c 

c 


7777S 


call  rore^r 

the  job  in  the  textport  needs  to  run  in  the  foreground 
call  KeepAsC  1024.767; 

the  window  neeas  to  have  an  aspect  ratio  of  a;b 
i  =  winope(  *Vorticity  Tool ‘.14) 
call  wintit( 'Vorticity  Tool'. 14) 
get  a  window  and  title  it 
call  getsizi  Xsize.  Ysize) 

find  out  how  large  it  is,  then  laake  the  cransf orcation 
matrix  which  will  scale  everything  appropriately. 

if(  Ysize. eq. 1024.0)  Ysize  =  767. 

mex  isn't  running.  Ysize  is  really  767... 
scale  -  1024.  /  Ysize 

1024.  IS  the  largest  X  or  Y  value  my  data  will  ever  nave. 


do  77777  1=1,4 
ao  77777  j=i ,4 

ifC  i.eq.j/  cnen 

i.hi.i;  = 

else 

UH  1 .  J  /  = 

endit 

cont  ir.ue 
m^4,4)  =  scale 


1.0 


0.0 


11^.2, 2;  =  -mv2,2) 
m^4,2)  =  Ysize 

this  shold  flip  the  image  vertically 


call  multtiaC  ;.u 
call  color!  -llllTE', 
call  clear 
uo  7777S  1=6,04 

call  -lapcoiC  i,  i*4,  J.  i50-',2*i)) 

cunt inue 


c 

c 

c 

c 

c 


ilere  are  the  preaefined  colors  1-7  (defined  in  'fgl.n') 


3  LACK  liED  GREEN  YELLOW 

BLUE  MAGQ)T  'CYAN  WHITE 


.i 
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□  define  UPSIDEJIOWN 
cO  define  MFRC 

c  <M>ichael  <F>iiia  <R>eader  <C>oords 

0  define  YYMIN  70.0 

c  actually  -YYMIN  ..  Y  =  Y  ^  YYMIN 

U  define  XXMIN  30.0 
c 

□  define  YSPLIT  5.217881548 

G  define  XSPLIT  5.217881548 
0  The  conversion  is: 

c  Y  =  CYofrc  f  YYMINj  *  YSPLIT  and  siailarly  for  X. 

subroutine  getfrmC  frame! >  frame2.  UuDim.  ierr) 
c 

c  Description: 

c  Read  in  the  data  files,  storing  the  points  lu  FRAMEl  and  FRAIiE2. 

c 

c  Return  Values. 

c  X.Y  points  in  frame!  and  frame2 

c  ierr  =  0  if  no  error  occurred 

c 

c  History: 

c  ^cpg>  aug  87  ---  vl.O 

c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

J  include  '"vorticity . h" 
c 

character*80  filnaml,  filnam2 
real  Xoffset.  Yoffset 
integer  NumPts 
external  NumPts 
logical  UnOim 
o 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c  INITIALIZE  both  framel  and  frame2 

c 

ao  1  1=1,  MaxPts 

frameK  i.X)  =  0.0 
framel 1  i.Y)  =  0.0 

frame2(  i.XJ  =  0.0 
frameZC  :.,Y)  =  0.0 

1  continue 

c 

c  FRAMEl 

1000  format! '  Please  enter  the  name  of  the  file  which  contains'/ 

T  '  the  data  for  frame  1.') 

1010  lormatC'  Error  opening  the  data  file  ' ,ao0/ '  Please  try  again.'/) 

10  conttnue 

c  return  here  on  error  Gl 

filnaml  =  'framel.dat ‘ 
write(  *,  1000} 
call  getline(  filnaml 1 
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c 

call  nalcefrinC  framel,  filnaml,  ierr) 
n  =  NumPcsC  framel •  2} 
call  CO  lor (  BLUE) 
do  15  i=l,  n 

call  circ(  framel(i.X),  framel(itY)t  3IGRAD) 

15  continue 

i£(  ierr.eq.l)  then 

c  cry  again  if  error  opening  the  frame  data  file 

writeC  1010}  filnaml 
goto  lU 

else  if(  ierr.ne.O;  then 

c  calling  routine  will  have  to  deal  with  the  error 

return 

end  if 

c 

c  fRAIlE2 

1020  tormatl'  Please  enter  the  name  of  the  file  wiiicn  contains'/ 

-  ■  tne  data  for  rrame  2.\ 

20  continue 

filnam2  =  'frame2.aat' 
write!  *,  1020) 
call  get  line!  filiiam2) 
c 

call  inattefrm!  frame2,  rilaam2,  ierr) 
if(  ierr.eq.l;  then 

write!  *,  1010)  £ilnam2 
goto  20 

else  if!  ierr.eq.O)  then 

e  get  the  X-  and  Y-o£fsets  for  this  frame 

write!  *.  1030) 

read!  *,  *)  Xoffset,  Yoffset 

i:  -  i.’umPtsv  framei.  2, 

uo  30  1=1,  n 

frame2!  i,X)  =  frame2!i,X/  -  Xoffset 
frame2!  i,Y)  =  fraTae2!i,Y)  -  Yoffset 

30  continue 

call  color!  CYAW/ 
ao  40  i=l,n 

call  circv  frame2!i,X),  frame2!i,Y),  BIGRAD) 
40  continue 

c 

call  convert!  1.0,1 .0,0 .0,0.0, UnDim) 
return 

end  if 

1030  format!'  Please  enter  the  X-  ana  Y-offsets  for  frame  2.'/ 

They  should  be  chosen  such  that  !0,0)  in  frame  1  is  the'/ 

point  (^Xoffset,  Yoffset)  in  frame  2.') 

end 

subroutine  makefrm!  frame,  filnam,  ierr) 
c 

c  open  the  frame  cata  file 

c  throw  away  the  first  (uescriptor;  tine 

c  wnile  not  EOF 

c  find  out  how  nany  points  are  in  the  next  set 
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c 


I 


c 


c 


read  X  values 
read  Y  values 

convert  X  and  Y  to  a  scd  coord  system 
end  while 


c 

c  The  following  is  an  example  of  the  data  contained  in  the  file: 


c  ORIGINAL  OBJECTS 

c 

c  OBJECT  LABEL 

0  1 

D  2 

D  3 

D  4 

D  5 

c  X-CENTROID 

225.75 

662.50 

497  .44 

861.73 

705.30 

c  Y-CENTROID 

109.25 

115.00 

133.78 

185.91 

234.30 

c 

c  OBJECT  LABEL 

Q  36 

0  37 

D  38 

c  X-CENTROID 

445.7  8 

570.09 

6  85.25 

c  Y-CENTROID 

916.00 

922.36 

940.25 

C123456739  123456789 

123456789 

123456789 

123456789 

123456789  123456789  12345 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

Li  include  "vorticity .h*" 

c 

integer  sett  start 
character*80  filnam 

character*80  linel.  Iine2.  Iine3.  junk 
integer  HowMany 
external  Uowl-Iany 
c 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

open(  3t  err  =  2000t  file=filnam«  status  =  'old') 
print*,  '  Reading  data  from  '•filnam 
c 

1000  formate  aSO) 

iOlO  formate  17 x,  5e4x,  £7.2 J) 

1020  formate'  Zero  objects  found  in  coordinate  set  ',i2,'.') 

read!  3,  1000,  end=2010}  junk 
c 

start=0 


100 

c 


c 


c 


start  stores  the  current  number  of  points  chat  have 
been  read  in. 

set  =  0 
c  ont inue 


reaae 

3. 

1000, 

fend=2020) 

juDK 

reade 

3. 

1000. 

end=2050) 

linel 

reade 

3. 

1000, 

end=2050) 

liae2 

read! 

3. 

1000. 

end=2050) 

line3 

set  = 

set 

;  +  1 

n  =  HovManye  linel) 
ife  n.gt.O)  then 

reade  line2,  1010)  eframe(i.X),  i=start-rl,  start*n) 
reade  line3,  1010)  eframeii.Y),  i=start+l,  start-^n) 
start  =  start  t  n 
else 

writee  *,  1020)  set 


CO 
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end  if 
c 

goto  100 
c 

2000  c  ont inue 

c  error  opening  the  frame  data  file 

ierr  =  1 
return 
c 

2010  continue 

c  empty  file 

ierr  =  2 
close(  3) 
return 
c 

2020  continue 

c  end  of  file  -  OK 

close(  3} 
ierr  =  0 

d  ifdef  HFRC 

ao  2025  x=l. start 

framed, Y)  =  (framed.Y,'  -r  YYMIU)  *  YSPLIT 
framed.X)  =  {frame(i,X)  +  XXMIN)  *  XSPLIT 

2023  continue 

end  if 

ifdef  UPSIDE_DOWK 
do  2030  1=1, start 

framed.Y)  =  -framed.Y)  +  1024. 

2030  continue 
ill  end  if 

return 

c 

2050  continue 

c  vrong  number  of  lines,  three  should  appear  at  once  -- 

c  OBJ  tIUMBEIR  /  X-COORDS  /  Y-COORDS - >  Frame  probably  incorrect 

ierr  =  3 
close(  3) 
return 
end 

integer  function  HowMany(  line) 

c 

c  Check  that  the  first  word  is  OBJECT, 

c  Then  count  to  see  how  many  u  signs  there  are. 

c 

c  We’re  parsing  this  line: 

c  OBJECT  UBEL  G  36  U  37  □  38 

C123456789  123456789  123456789  123456789  123456789  123456789  1234567S9  12345 
c 

c  so  we'll  only  look  every  11  places  for  the  G. 

c 

characterise  line 
integer  pos 
characteri6  word,  obj 
character*!  lb 


c 
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Do 
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integer  function  NumPtsC  points*  dim) 
c 

Li  include  "vorticity .h"’ 
c 

integer  dim 

real  points(  MaxPts*  dim) 
c 

n  =  0 

10  continue 

c 

if(  points(n-*-l .  X}  .eq.0 .0  .and.points(n-t-l  >  Y).eq*0.0)  goto  20 
n  =  n-*-! 

i£(  n.eq. MaxPts)  goto  20 
goto  10 

c 

20  continue 

c 

c  "n"  now  contains  che  number  of  points 

c 

MumPts  =  n 
return 
enn 

integer  function  iiumPoiyC  poLyj 

include  "vorticicy. h" 
c 

n  =  0 

100  continue 

if!  poly(n+i ,1 ; .eq.O )  goto  200 
n=n+l 

if(  n .eq.MaxPoly )  goto  200 
oOto  ICO 
c 

200  continue 

c 

c  "n"  now  contains  che  number  of  polygons 

c 

Ii'umPoly  =  n 

return 

end 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


subroutine  sortfrmC  frame) 

Description:  Uses  a  simple  exchange  sort  to  order  FRAME 
on  decreasing  Y.  Ignores  (OtO)  pairs. 

This  sort  should  be  efficient  enough  given  the  small  number 
of  points  which  can  be  in  a  frame.  If  MaxPts  goes  above  200 • 
someone  should  consider  exchanging  this  for  a  shell  sort  or 
something  else  more  efficient. 

Author:  Chuck  Gendrich 

History: 

August  •  1987  vl  .0 
Procedure: 

Determine  how  many  points  there  are 
For  all  points  except  the  last 

check  this  point  against  all  subsequent 
if  X  >  Xsubsequent. 

swap  X*s  and  Y's 
if  X  ==  Xsubsequent 

check  Y's  and  swap  if  appropriate 

end  for 


ccccccccclccccccccc2ccccccccc2ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc 

c 

i_i  include  "vorticity.h" 
c 

n  =  NumPts(  frame.  2) 

c  find  out  how  many  points  there  are  in  this  frame 

c 

uo  100  i=l,n-l 

c  check  all  points  but  the  last 

do  75  j=i+l.  n 

c  compare  it  against  all  remaining  points 

if(  f rame( i.Y) . It .frame( j «Y) }  then 

call  suap(  framed. X).  frame(  j.X)) 
call  swap(  framed. Y).  frame{  j.Y)) 
else  if  (frame( i.Y) .eq. frame! j ,Y) )  then 
c  compare  X's 

if(  frame(iiX). It. frame! j.X))  then 

call  swap!  frame! i.X).  frame!  j.X)) 
call  swap!  frame!!. Y).  frame!  j.Y)) 

end  if 

end  if 

75  continue 

100  continue 
return 
end 

suoroutine  swap!  a.  b) 


real  a.  b.  tmp 
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subroutine  HatchesC  framel.  frame2.  match,  ierr) 
c 

c  Description:  Establishes  a  correspondence  between  points 

c  in  framel  and  frame2.  Match  is  setup  such  that 

c 

c  framel (  i)  < — >  frame2(  match(i)) 

c 

c  Algorithm:  (What  algorithm?  Huh?  Who?  Where?  :-(  ) 

c  We'll  be  looking  for  points  within  a  certain  neighborhood  of 

c  the  original  point.  At  the  moment,  we'll  assume  there's  only 

c  one  point  which  will  fall  into  the  neighborhood,  but  that  may 

c  not  be  the  case, 

c 

c  We'll  error-check.  If  more  than  one  point  matches,  this 

c  routine  will  return  an  error.  If  no  point  matches,  a  zero 

c  will  be  put  in  the  appropriate  place  in  Match, 

c 

c  Author:  ChucK  Gendrich  (**I  wish  I  knew  how  to  do  this  one  better, 

c 

c  History: 

c  August.  1987  vl .0 

c 

ccccccccclccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc 

include  "vorticity .h** 
c 

real  close,  clo8e2 
paramecer(  close  =23) 
parameter(  close2  =  close^close) 

c  half  the  distance  in  pixels  between  the  two  closest 

c  points  in  the  original  image  (FENCED 

c 

real  dist 
integer  NumPts 
external  dist.  NumPts 
c 

integer  found 

c  how  many  matches  have  we  found 

c 

c  initialize  Match 

c 

do  10  i=l,  MaxPts 

match(i}  =  0 

10  continue 

c 

c  OK,  we'll  do  it  the  hard  way.  Compare  every  point  in  framel 

c  with  every  point  in  £rame2.  Tally  all  the  distances  less  than 

c  CLOSE.  If  there’s  only  one  point,  wonderful.  Otherwise,  oops, 

c 

nl  =  NumPts C  framel,  2) 
n2  =  NumPts (  £rame2,  2) 
c 

do  100  i=l,  nl 
found  =  0 
xl  =  frameKi.Xj 
yl  =  framel (i.Y) 
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75 


c 

c 

c 


100 


do  75  j=lt  a2 

x2  =  £rame2(j.X} 
y2  =  £rame2(j.Y} 
dx  =  x2-xl 
dy  =  y2-yl 

if(  (dx*dx  +  dy*dy) .lC.clo6e2)  then 
£ound  =  £ouad-t-l 
i£(  found. gC.MacchOiis)  then 
ierr  =  6 
return 

end  if 

matchs(  found)  =  j 

end  if 
continue 

i£<  found. eq.O)  then 
macch(i)  =  0 

else  i£(  found.eq.l)  then 
matchCi)  =  matchsC  1) 

else 

Of  course  this  should  oe  made  more  sophisticated  so  that 
it  can  conflict-resolve  any  problems.  But  for  now»  Just 
chalk  up  an  error  and  return, 
ierr  =  7 
return 


end  if 
continue 
ierr  =  0 
return 
end 
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subroutine  VelocitiesC  frameli  £rame2.  match,  vel) 
c 

c  Description:  vel_mker  computes  the  position  and  botn 

c  components  of  the  velocity  vector  between  matching  points 

c  in  framel  and  f rame2 . 

c 

c  The  position  is  assumed  to  be  at  the  midpoint  between  the  two 

c  matching  points.  If  there  is  no  matching  point  in  FRAME2  for 

c  some  point  in  FRAliEl  (.match(i)  ==  0),  the  position  and  velocity 

c  components  are  left  0. 

c 

c  The  initial  "velocity"'  is  in  dPixels.  A  conversion  for  pixels 

c  to  real  units  of  measure,  and  knowledge  of  dTime  between  frames 

c  is  required  so  that  we  can  calculate  the  real  velocity.  This 

c  will  be  uone  in  a  subsequent  routine, 

c 

c  Autnor;  Chuck  Gendrich 

c 

c  History; 

c  August.  19  87  vl.O 

c 

ccccccccclccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7c 

c 

u  include  “vorticity .h" 
c 

characcer*80  line 

integer  HumPts 

real  avg,  getreal 

external  avg.  IlumPts.  getreal 

c 

c  initialize  the  velocity  descriptor  array 

c 

uo  10  1=1.  MaxPts 

vel(i.X)  =  0.0 
vel(i.Y)  =  0.0 
vel(i.Vx)  =  0.0 
veKi.VyJ  =  0.0 

10  continue 

c 

writet  *.  1000) 

1000  format!'  Please  enter  now  many  pixels  a  point  moving  at’/ 
tne  free  stream  velocity  will  move  between  frames.'/ 

!Vx)  =  CdeltaX)  +  baseV) 
line  =  '  -36.1655  pixeis/f rame' 

baseV  =  getreal!  line; 
scale  =  4.0  *  exp!  -absibaseV)  /  28.0) 
c 

a  =  N’umPts!  framel.  2) 

1  =  0 
c 

do  100  j=l .  n 

c  laake  a  corresponding  entry  in  vel  for  eacn  point  in  framel 

if!  match! j ; .ne .0 )  then 

c  we  have  something  to  work  v/ith 
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i  =  i+1 

vei(i,X)  =  avg(framei(  j.X),  frane2(  raatchC  j  /  .X) ) 
veKi.Y)  =  avg(framel(  j.Y),  trame2(  matchtjj.Y)) 
vel(i.Vx)  =  frame2(niatch(  j)  .X)  -  trameKj.X)  ■*-  baseV 
velCi.Vy)  =  traine2(macch(  j)  ,Y)  -  irameKj.Y) 
call  laovali.  vel(  i,X),  vel(i,Y)) 
call  vectorV  vel(  t.Vy),  veK  i,Vx)i  scale) 
c 

endif 

100  continue 

c 

c  Now  construct  the  TURB3d  output  files, 

c 

opent  unit=3.  rile= 'vel.xyz' ) 

c  unit3  —  flow  field  point  locations 

open(  unit=4t  file= 'vei .q ' ) 

c  unit4  —  tlow  tiela  information 

c 

write!  3.  lOlU;  i 

c  write  out  tne  X  and  Y  locations 

write(  3.  1020)  vveK  j.X),j=l,i/ 

write(  3t  1020;  ^vei(  j,Y),j=l,i/ 

c 

write!  4,  1030)  i 

c  'write  out  the  Vx  and  Vy  velocities 

write(  4,  1020)  (veK  j,Vx),j=l.i) 

write(  4,  1020)  >vel(  j.Vy),j=l,i; 

do  110  j=l . i 

c  write  cut  a  fictitous  pressure  value 

write!  4,  1040) 

110  continue 

do  120  j=l,i 

c  write  out  a  fictitous  vorticity  value  tfor  now; 

writer  4,  1040/ 

120  continue 

1010  format!  i3,'  1’) 

1020  iormatC  5(fl0.5ilx/) 

1030  formate  i3 , '  I'/'l  1  1  1') 

1040  rormat!'  1') 

c 

return 

ena 
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real  function  distC  xl.  yl .  x2.  y2) 
c 

ccccccccclccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc 

c 

real  xl.  yl ,  x2,  y2.  sqrd 
c 

sqrd  =  (xi-y-2)**2  (yi-y2}**2 

c 

disc  =  sqrt(  sqrd} 
c 

return 

end 

real  function  avg(  a>  b) 
c 

ccccccccciccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc 

0 

real  a.  b 

c 

avg  =  Ca+Dj/2.0 


return 

eno 


real  function  tanl(  ay,  dx) 
c 

ccccccccclccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc 

c 

c  returns  the  arccan  in  ccgrees  given  dx  and  dy 

c 

real  dx.  dy,  angle,  pilCO,  small 
c 

parameter!  small  =  0.0001} 

0  to  xeep  rrom  diviaing  oy  aero 


c 

c 


parameter!  pilSO  =  57.295779513} 

180/pi  ror  converting  to  degrees 

if!  abs!  dx} . It.smaii}  then 

let's  not  aivide  by  zero,  shall  we? 
if(  dy.gt.O}  then 

tanl  =  90  .0 

else  if(  dy.lt.U)  then 
tanl  =  270.0 

else 

tanl  =  0.0 
???  punt  ??? 

endif 

else 

angle  =  atani  ay/ox}  *  prloO 
ifC  dx.gt.O)  then 

rirst  or  fourth  quadrant  --  atan  is  OK 
angle .  It  .0 .0 /  then 

angle  =  angle  +  360.0 
make  sure  angle  is  >  0 


enuif 


Dec  9  23  ;  52  1  987  ^^eometry.f  Page  2 


else 

tanl  = 

angle 

c 

endif 

second 
tanl  = 

or  third  quadrant  --  add  180  deg. 
angle  +  180.0 

eadif 

return 

end 


real  function  interp(  xxO.  ffxO.  xx2.  ffx2»  xl.  ierr; 
c 

ccccccccciccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc 

c 

c  returns  the  interpolated  value  for  fix)  given  the  function 

c  value  at  boundaries  xO  and  x2.  xl  has  to  lie  between  xO  and  xl . 

c 

real  xO,  fxO>  x2<  fx2«  xl,  xxO.  tfxO,  xx2.  i;£x2 
integer  ierr 

c 

real  soall 

parameterv  scail  =  0.00001; 
c 

real  dxO,  d£x 
c 

xO  =  xxO 

fxO  =  ffxO 

x2  =  xx2 

£x2  =  £fx2 

c 

c  since  we  uon’t  want  to  change  the  original  values  for  xO  and  x2, 

c  etc.«  we  tia^e  our  own  copies  wnich  we  can  switch  around  if 

c  necessary, 

c 

inter p  =  rxO 


c 

c 

c 


t 


c 


c 

c 

c 


(.junk  answer  just  in  case  they  don't  check  ierr  and  u'.terp's 
value  isn't  changed  because  of  some  error... 7 

dxO  =  x2  -  xO 

ifi  ansi  dxOJ  .it.smalU  then 

let's  avoid  division  by  zero,  eh? 

ierr  =  16 

return 

else  if(  dxO.lt.O)  then 

xO  >  x2  - >  Ve  have  to  switch  them  around 

call  swapC  xO.  x2) 
call  swap!  fxO,  fxZi 
dxO  =  -dxO 

else 

they're  already  in  the  rigat  order  anu  plenty  lar  apart 

enaif 

if(  xi  .gt.x2.or.:jl.lt.xO;  then 
out  01  bounds  (oops; 
ierr  =  17 


return 
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end  if 
c 

dxl  =  xl  -  xO 
dfx  =  fx2  -  fxO 
c 

interp  =  fxO  +  (  dfx  *  dxl  /  dxO) 

return 

end 

subroutine  centroidC  poly,  vel,  i.  xO,  yO,  ierr) 
c 

c  Description:  The  centroid  of  a  four-sided  figure  is 

c  calculated.  This  is  located  at  the  intersection  of  the  two 

c  lines  which  join  the  midpoints  of  opposing  sides, 

c 

c  Author:  Chuck  Genorich 

c 

c  History: 

c  August,  1987  vl.O 

c 

ccccccccclccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc 

c 

cLi  uetine  DEBUG 
u  include  "vorticity.n” 
c 

real  small 

paraneterC  small  =  0.00001} 
c 

real  avg 
external  avg 


c 

c 

c 

c 

c 

c 

c 


real  xO,  yO 

the  coordinates  for  the  centroio  A 


we'll  define  the  polygon  as  one  consisting  of  \ 
four  corners.  A,  3,  C,  and  D.  The  coordinates  \ 
of  tne  corners  are  Ax,  Ay,  Bx,  By,  etc....  a,b,c,  \ 
ana  d  point  to  the  appropriate  entries  into  vel.  D. 


integer  a,b,c,d 

real  Ax,  Ay,  3x,  By,  Cx.  Cy,  Dx.  Dy 


3 


real  ml,  m2,  ol,  b2 

c  we'll  calculate  the  slope  and  intercept  of  both  lines 

c  which  join  the  midpoint  of  opposing  sides.  To  find  the 

c  intersections,  we  let  yl  =  y2  and  solve  for  x.  Doing 

c  this  we  find  xO  =  -(b2  -  bl}/(m2  -  ml}.  Plugging  oack 

c  in ,  we  f  ind  yO  =  ml  *  x  +  b  1  t  =  nL2  *  ;c0  +  b  2 . . . . } 

c 

real  deitaY,  deitaX,  yl ,  xl,  y2,  x2 
c 

ccccccccclccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc 

c 

a  =  poly(  i,ly 
b  =  polyl  i,2} 
c  =  poiyi  1,3} 
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d  =  poly(  i,4) 
c 

Ax  =  vel(  a.X) 

Ay  =  vel(  a,Y) 

Bx  =  vel(  b,X) 

By  =  vel(  b,Y) 

Cx  =  vel(  c.X) 

Cy  =  vel(  c,Y) 

Dx  =  vei(  d,X) 

Dy  =  vel(  d.Y) 
c 

u  ifdet  DEBUG 

wrice(  *,  990)  a>Ax>Ay>b«Bx.ByiC.Cx.Cy,d.Dx>Dy 
990  formate/'  centroia;  finding  the  centroid  with  these  corners 

-  (4x.i2,':  ('.f8.3.'.',f8.3.')')) 

1000  format('  V/ierd  centroid  U.o  cx  or  dy)  at  ’  ,f  3.3  , ' .  ’ .  £8 .3  , ' ) ' ) 

1010  formatC*  centroid;  neither  line  vertical.  nl :  ',£8.3, 
m2;  '.£8.3/'  bl;  '.£8.2.'  b2;  ',£8.2.'  — >  ('.£8.3.'.'. 

-  t8.3.')'/) 

1020  format('  centroid:  ’,a6,'  vertical,  m:  ',£3.3,'  b:  ', 

-  £8.3/’  -->  ('.£8.3. '.'.£8. 3,')'/) 

i_  endii 

c 

c  —  initialize  xO  ano  ierr  (xO  must  oe  known  later; 
c 

xO  =  -99999.0 
lerr  =  0 
c 

c  —  line  1  first 
c 

yl  =  avg(  Ay,  By) 
xl  =  avg(  Ax,  BxJ 

v-elcaY  =  avgt  Dy,  Cyw  -  yl 
deltaX  =  avg(  Dx,  Cx)  -  xl 
c 

i£(  abs(  deltaX) . It. small)  then 
c  let's  see  i£  we'll  be  dividing  by  zero 

if(  absl  ueltaY) . it. small)  then 

c  inidpt(AB)  =  miopt(CD)  i?)  and  the  centroid  is  there 

yO  =  yl 
xO  =  xl 
ifdet  DEBUG 

write!,  *,  1000)  xQ,  yO 

endif 

return 

else  i£(  abs(  deltaX/deitaY ), It . small)  then 

raidpt(AB) - >  midptiCD;  is  a  vertical  line 

mi  =  xO 
xO  =  xl 

end  if 

else 

ml  =  deitaY  /  deltaX 
ol  =  yl  -  iai*xl 

enaif 

c 
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c 

c 


c 


c 


c 


u 

u 


c 

c 

c 


c 

c 

c 

c 


L. 

b 


c 


line  2  next 


y2  =  avg(  Ay.  Dy) 
x2  =  avg(  Ax.  Ox) 

deltaY  =  avg(  By.  Cy)  -  y2 
deltaX  =  avg(  Bx.  Cx)  -  x2 


i£(  abs(  deltaX) .It. small)  then 

i£(  abst.  deltaY)  .It. small/  then 

midptCBC;  =  midptCDA)  (?)  and  the  centroid  is  there 
yO  =  y2 
xO  =  x2 


else 


encif 


ifdef  DEBUG 

write(  *,  1000)  xO,  yO 

endif 

return 

else  i£(  abs(  ueitaX/ueitaY/ . It . smal 1)  then 

aidpt(BC)  - >  cidptlDA)  is  a  vertical  itne 

ifv  xO .ue . -99999 .0 ;  then 

Che  ocher  line  was  also  vertical 
Punt ! 
ierr  =  4 
return 

else 

m2  =  xO 
xO  =  x2 

end  if 

enoif 

ci2  =  ueltaY  /  deltaX 
b2  =  y2  -  ra2*x2 


aow  tor  the  intersection  vi.e.,  the  centroid) 

calculate  xO  i£  we  still  need  to 
ift  xO.eq. -99999.0)  then 

neither  line  was  vertical 
xO  =  -(62  -  bl)  /  {m2  -  tall 
yO  =  a2  *  xO  +  b2 
ifdei  DEBUG 

write(  *,  1010)  ml,  ia2,  bl,  b2,  xO.  yO 

endif 

return 

endif 


calculate  yO  using  the  appropriate  slope  ana  intercept. 
if(  mi .eq. -99999.0)  then 

line  1  is  vertical  so  use  u2  anu  b2 
ifdef  DEBUG 

writeC  *,  1020)  'line  1',  ia2,  d2,  xO,  yO 
endif 

yO  =  u2  *  xO  +  b2 

else 

use  ml  ana  bl 
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subroutine  Folygonsv  vel.  poly,  ierr) 

Description:  Establishes  a  polygon  descriptor  array.  Each 

element  in  polyCi)  describes  a  corner  of  the  polygon  according 
to  the  following  scheme: 


poly(  i.l) 
polyv  i>2j 
poly(  i,3) 
poly(  i.4) 


— >  upper  LH  corner 
— >  upper  RK  corner 
-->  lower  LH  corner 
— >  lower  RH  corner 


c  Procedure:  Try  to  make  roughly  rhomboidal-shaped  polygons.  Try 

c  to  use  each  point  in  VEL  as  an  upper  LH  corner.  If  it  doesn't 

c  work,  move  on  to  tde  next  point, 

c 

c  Author:  Chuck  Genurich  ^"SCORPIO;  ail  Scorpio's  are  assassinateo 

C 

c  History: 

c  15  aug  87  <cpg>  vl.l  —  removeu  the  triangle  logic.  look  in 

c  oiaer  versions  to  recover  it. 

c 

ccccccccclccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc 

c 

L)  include  “vorticity.h" 
c 

real  delta 

c  when  we  look  for  a  corresponding  corner,  it  will  be  in  a 

c  particular  direction  +/-  delta  deg. 

c 

integer  ULH,  URH,  LRH,  LLH 

c  pointers  to  the  UpperLef tUand.  UpperRK,  lilH,  and  LLH 

c  vertices  of  the  polygon  under  construction 

c 

logical  IncrP 
integer  KumPts,  look 
external  ITumPts,  Look,  IncrP 


n  =  NumPts(  vel.  4; 
delta  =  30.0 

look  for  point  C  direction  +/-  delta  aeg. 

p  =  0 

we  have  currently  found  zero  polygons 
do  100  ULH=1,  n 

URH  =  look(  ULH,  vel,  4,  RIGHT,  delta,  imtchs.  lerr; 

if(  URH. ne. BAD. and . ierr .eq.O)  then 
found  tne  URH  corner 

LRH  =  look(  URH,  vel,  4,  DOWN,  delta,  watchs,  ierr; 
if(  LRH.ne. BAD. and. ierr. eq.O)  then 

LLH  =  look(  LRK,  vel,  4,  LEFT,  delta.  ..atcns,  ierr; 
ifl  LLH.ne.BAD.ana. ierr .eq.O)  then 
we  have  an  OK  polygon 
if(  IncrP I  p;)  then 
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c 


c 

c 


c 

100 


poly(  p.  1)  =  ULH 

poly(  p.  2)  =  URH 

poly(  p,  3)  =  LRH 

poly(  p.  4)  =  LLH 

else 

return 

doesn't  o^ke  sense  to  look  for  any  uiore 
polygons  if  poly  is  full.... 

endif 

end  if 
enaif 
endif 

if(  ierr.ne.O/  then 

there  was  an  error  "looking" 
return 

endif 

in  every  ocher  case,  we  couldn't  finn  a  sacisiactory  polygon 
Forget  using  tnis  corner  as  the  start  of  a  polygon. 

cont inue 

return 

end 


logical  function  IncrP(  p} 
c 

c  s imple-mindeo  routine  to  increment  p  or  quit  if  there's  a 

c  problem 

c 

c  returns:  TRUE  if  p  is  <=  tne  size  of  poly 

c 

ij  inciuae  vorticity.n" 
c 

p  =  p-^l 

IncrP  =  ip .  le .liaxPoly ; 


c 


return 
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c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


integer  function  look(  p. points. din. angle»delta,  matchs,  ierr) 

Description:  "look”  uses  the  location  of  point  '?'  among 

the  given  points.  It  examines  all  the  other  points  tor  ones 
which  are  between  100  and  195  pixels  away,  lying  between 
vectors  at  angles  "angle  +  delta"  and  angle  -  delta  . 


->>  angle  +  delta 


_-=/  no  good 

_-=/  ^  ) 

_-=/  )  ) 
p  -  100  195 

o  < - I - I - >  angle 

)  ) 

\-=_  )  OK  J 

\-=_  >  .  / 


\-=_ 

\->>  angle  -  delta 

Return  value:  "Look"  returns  SAD  if  no  point  tails  within 
the  defined  search  region. 

If  only  one  point  matches,  "look"  returns  the  pointer  to  that 
point . 

If  niore  than  one  point  matches,  "look  returns  a  pointer  to  the 
one  wnich  is  closer  to  the  preferred  direction. 

Autnor:  Chuck  Genorich 

History : 

August.  1987  vl  .0 


ccccccccclccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc 

c 

^  include  "vorticity . n" 

c 

integer  dim 

real  pointsl  liaxFts.  aim; 


real  angle,  delta 

integer  NumPts 
external  NumPts 


I 


c 


integer  founc 
real  ApIusD.  AminusD 

--  to  score  /uiglo  Delta  anu  Angie  -  Delta 

real  ::1.  jl.  x2.  ^1.  a;.,  ay 
real  d.  direction 

distance  and  airection  from  vxl.yi;  to  1x2, y21 


real  I-linD.  MaxD 
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parameCerC  MinD  =  40.0) 

parameter (  MaxD  =  180.0) 

c  the  closest  and  farthest  that  a  point  has  to  be 

c  in  order  for  it  to  be  considered, 

c 

real  Distnc(  MatchDim) 

c  difference  between  Che  desired  and  the  actual  angle 

c 

integer  best 

c  a  pointer  to  the  best  smtch  (when  there's  niore  than  one) 

c 

logical  first.  DoDraw 
character*!  char 
save 

data  first/ .TRUE./ 
c 

ccccccccciccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc 


if(  FIRST)  tnen 

first  =  .FALSE. 

print*. 'Draw  in  ’"look**  lines?  [n]* 
read( ' (Al ) ' )  char 

DoDraw  =  (char. ne. 'n' .and. char. ne. 'N' .and. cnar.ue. '  ') 

end  if 


lerr  =  0 

n  =  NuaPcs(  points,  aim) 
do  1  1=1.  MatchDim 

matchs(i)  =  0 

cont inue 

if(  n.lt.p)  then 

i.e.,  are  we  oOing  to  be  doing  something  reasonaole? 
tf  we  ^et  here,  the  answer  is  "iiO” 
terr  =  11 


return 


end  if 


round  =  0 

AplusO  =  angle  aelta 
AminusD  =  angle  -  delta 

Get  tne  X  and  Y  coords  we're  going  to  be  comparing  against,  then 
look  at  every  otner  point... 

xl  =  points!  p.  X) 
yl  =  points!  p.  Y) 
if(  DoDraw)  then 

call  CO  lor (  WHITE) 

call  circ!  xi.  yl ,  JiGf;AD) 

call  circ(  xl.  yl ,  SMALLRAD) 

encif 


CO  100  1=1.  n 

...to  avoid  excessive  indenting  and  nesting,  tne  tests 
will  be  for  disqualification,  if  a  point  disqualifies 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


100 

c 

c 

c 


c 


ve’II  simply  'concinue'  on  Co  Che  nexc  one. 

i£(  i.eq.p)  goco  100 

x2  =  pointsC  i,  X) 
y2  =  poincs(  i,  Y) 
i£(  DoDraw)  Chen 

call  circ(  x2>  y2.  SMALLRAD) 

endi£ 

dx  =  x2-xl 
dy  =  y2-yl 

d  =  sqrc(  dx*dx  +  ciy*dy) 

i£(  d. IC.MinD.or .d.gt.MaxD)  goCo  100 

direction  =  canl(  uy,  dx) 

if(  utrecCion.oC.AplusD)  direccior.  =  cireccion  -  360.0 
MB  --  tanl  always  returns  0  <=  angle  <=  360 
but  i£  AminusD  <  0  and  the  direction  really  *ts* 
OK  at  its  position  in  the  fourth  quadrant,  it'll 
never  i.atch.  Plus,  we  have  to  have  A+D  >  A  -  D... 
Therefore,  i£  direction  >  A+D,  we'll  give  it  a 
chance  to  match  with  angles  in  the  Ath  quad. 

i£(  direction. le.AplusD. and. direction. ge. AminusD)  Chen 
found  =  found  +  1 
i£(  found. gc.MaCchDim)  then 
ierr  =  12 
return 
enoif 

if(  DoDraw)  then 

call  circv  x2,  y2,  BIGDAD) 
call  colorC  RED) 
call  raove21  xl,  yl ) 
call  draw2(  x2,  y2) 
call  color!  ’.JHITE) 
end  if 

Distnc!  found)  =  dist(  xl,  yl ,  x2,  y2) 
macchs!  lounu)  =  i 
end  if 

cone inue 

So  now  v/e've  looked  at  every  point 

if(  founa.eq.O)  then 
look  =  BAD 

else  if(  found. eq.l)  then 
look  =  iiiatcns!!; 

else 

look  for  the  match  which  deviates  the  least  from  the 
desired  angle 
best  =  1 

do  200  1=2,  found 

ifC  DiseneCi) .Ic.Discnclbesc))  then 
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subroutine  fluids(  vel,  poly,  UnDim,  ierr) 
c 

c  Description:  calculates  the  fluid  mechanical  properties  of  the 

c  flow.  The  object  is  to  produce  pLot(s)  of 

c  (1 )  du/ dy 

c  (2)  du/dx 

c  ( 3 )  dv/ dy 

c  (4)  dv/dx 

c  ( 5 )  uv 

c  (6)  Wz  (vorticityj 

c 

c  Output:  ail  output  will  oe  v;ricten  to  stdout.  If  you  want  it  in 

c  a  file,  put  a  tee  on  the  process  when  you  run  it. 

c  e.g.:  freyja>  vorticity  1  tee  output_file 

c 

c  Procedure:  First  the  conversion  factor  for  pixel  <-->  mm  is 

c  obtained,  along  with  the  time  between  each  frame  and  the 

c  absolute  Y  offset  for  ooth  frames, 

c 

c  Author:  Chucx  Gendrich 

c 

c  history: 

c  August,  1987  vl .0  —  lots  of  linear  interpolation,  ugh. 

c  15  aug  87  <cpg>  vl.l  --  removed  triangle  logic,  see  older 

c  versions  to  recover  it. 

c  16  aug  87  <cpg>  vl.2  —  added  the  non~dim  title  stuff 

c  22  jan  88  <cpg>  vl.3  —  graphics  to  draw  the  polygons 

c  7  apr  88  <cpg>  vl.4  --  defines  a  colormap  based  on  h’zCnin-max] 

c 

ccccccccclccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc 

c 

cLi  derine  DEBUG 
^  include  "vorticity.n" 

real  ud,  d,  avg,  aist 

integer  HumPts,  nuraPoly 

external  NumPts,  b’umPoly,  u.  avg,  uist 

integer  a,  b 

c  a  ana  b  point  to  the  vertices  of  the  side  along  wnich 

c  the  velocity  is  to  be  integrated  next 


real  suba,  subi 

real  area,  intgrl.  GAMliA 

c  area  of  the  polygon  whicn  is  integrated 

c  integral  of  V  *  ds  arouna  the  polygon 

c  GAMMA  IS  the  total  circulation  tnroughout  the  frame 


c 


c 


logical  UnDim 

iroin  convert...  lEUE  if  results  are  oeino  UnDim'a 

cnaracter*2  .lum 
cnaracter*5  Afunits,  result 
cnaracter*/  Runits 

for  titling  the  output  vuepending  on  wnether 

the  results  are  oeing  uon-aimensionalizea  or  not... 
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real  Parray(  2,  4) 

parray  ■  points  array>  XY  locations  of  polygon  vertices 
iarray  -  intensity  array.,  vorticity  at  the  corners. 

real  Wzs(  HaxPolylt  VZmax,  tCmin 

Stores  the  values  of  Viz.  their  i.iax>  and  min. 
integer  Csplit,  shade 
parameter (  Csplit=  56) 

the  total  number  cf  colors  which  we  can  draw 


f ormat( ' 


X'  .a5,' 


Y  *  ,a5 ,3x,a5 ,a7 ) 


p  =  NumPoly(.  poly) 
n  =  I5umPts(  vel.  4) 
if(  UnDira;  then 

XYunits  =  '(noD;' 

else 

XYunits  =  'immJ  ' 

end  if 

:==  vorticity 

if(  GoOnO)  tnen 

title  =  'vorticity' 

end  if 

result  =  ■  Wz  ' 
call  CutHere^  title) 
if(  UnDim)  tnen 

Run ICS  =  '  ( noD  / ' 

else 

Runits  =  '(1/ sec ) ' 

end  if 

write(  *,  1000)  ATunics,  XYunits .  result,  Runits 

GAUIIA  =  0.0 
call  color(  UllITE, 
call  clearO 
WZmin=+999999. 
wZinax=-999999. 
do  100  i=l,  p 

call  color(  RED) 
area  =  0.0 
intgrl  =  0.0 

do  50  j=l,3 

a  =  poly(i,j) 
b  =  poly  (  i.  j-^1 ) 

^  and  b  are  tiie  enapoints  of  the  side  aion^  waicii 
V  *  ds  is  to  be  integrateo. 
call  integrtt  a,  b,  vel,  suba,  suoi) 

xi  =  veU  a,  a) 

X  2  =  V  e  i  V  u ,  A ) 
yl  =  veit  a,  Y) 


9 
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y2  =  veil  b.  Y) 
c 

call  ii:ove2(  xl.  yl ) 
call  draw2(  x2,  y2} 
area  =  area  suba 
intgrl  =  intgrl  +  subi 
c 

30  concinue 


c  now  for  the  last  side  (4  — >  1) 

a  =  polyli,4) 
b  =  poly(i.i) 

call  integrc(  a.  b»  vel.  suba.  subi) 

xl  =  vel(  a.  X) 

x2  =  vel(  b,  X) 

yl  =  vel(  a.  Y) 

y2  =  veil  b,  Y) 

c 

call  aiove2(  xl.  yl ) 
call  draw2(  >:2,  yl) 
c 

area  =  area  ^  suba 
intgrl  =  intgrl  +  subi 
Viz  =  intgrl  /  area 
Wzs(i)  =  Mz 

if(  Wz.  Ic.llZniia)  WZmin=Uz 
i£(  Wz.gt .VlZtaax)  WZu^x=Vlz 
c 

call  centroid(  poly,  vel,  i,  xO,  yO,  ierr) 
rf(  ierr.ne.O)  goto  100 
c 

intgrl  =  intgrl  *  50.0  *  .029520809  *  .029520809  /  1.0 
GAMMA  =  GAMMA  +  intgrl 
_  if net  DEBUG 

print*,'  Area ; ' ,area, '  intgrl; ' .intgrl, '  GAMMA  now. '.GAMMA 

u  enaif 

call  convert(  ;i0,  yO,  Wz,  0.0,  UnDim) 

call  CO  lor (  CYAN! 

call  cinov2(  xO,  yOl 

write(  num,'(i2}')  i 

call  charstC  num.  2) 


c 

100 

c 


c 

Lj  ifdei 
i_:  eiioif 


cont inue 
GAMMA  =  GAMMA 

circ  =  (pix**2  /  frame)  *  (frame/sec)  *  (mm/pix)**2 
print*. 'Total  circulation  is  '.GAMMA,'  (nun~2/sec ) . ' 
if(  GoOnO)  then 
do  150  i=l.p 

draw  shaded  polygons 

shade  =  vWzsii;  -  V.'Zain//(WZiwx-WZciia/*(Cspiit 

DEBUG 

print*, 'Vorticity  ',i,'  is  ’.V/zsli),'. 

call  colorC  shade! 
do  125  j-1,4 

ParrayV  X,  j;  =  vei(  poiy(i,j),X; 


Shade  is  '.shad 
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Parray(  Y,  j)  =  veil  poly(i,j),Y) 
continue 

call  polf2(  4,  Parray) 
cont inue 
end  if 

if(  GoOnO)  then 
continue 

end  if 
du/dy 

title  =  'du/dy' 
result  =  'du/dy' 
call  CutHereC  title) 

wrice(  *,  1000)  XYunits,  XYunics,  result,  hunits 
do  200  i=l ,  n 

dd  =  u(  Vx,  udy,  i,  vei,  tert , 

ignore  dd  ana  ierr  --  etcher  souething  was  rounc  c-nd 
printed  or  it  wasn't...  dd  produces  its  own  output. 

c  one inue 
au/dx 

title  =  'du/dx' 
result  =  'du/dx' 
call  CutHere(  title) 

write(  *.  1000)  XYunits,  XYunits,  result,  Runits 
do  300  1=1,  a 

dd  =  d(  Vx,  ddx,  i,  vel,  lerr) 

continue 
av/ dy 

title  =  'dv/ay' 
result  =  'dv/dy' 
call  CutHereC  title) 

writeC  *,  1000)  XYunits,  XYunics,  result,  Runits 
ao  400  i=l,  a 

dd  =  d(  Vy,  ddy,  i,  vei,  lerr) 

continue 
uv/  d:: 

title  =  'dv/dx' 
result  =  'dv/dx* 
call  CutHereC  ticle> 

writeC  1000)  XYunits,  XYunits,  result,  Runits 
ao  500  1=1,  n 

dd  =  aC  Vy,  udx,  i,  vel,  lerr) 


300 


cent inue 
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c  ===  Reynold's  stress 

c 

title  =  'instantaneous  Reynolds  stress' 
i£(  UnDim)  then 

result  =  '  uv  ' 

else 

result  =  'uv  ^m' 

Runits  =  'm/s)”2  ' 

end  if 

call  CutUere(  title) 

write(  *,  1000)  XYunits,  XYunits.  result.  Runits 
c 

do  600  i=l.n 

call  convert(  veHi.X),  vel(i.Y),  0.  velC  i,  Vx)*vel(  i ,  Vy) , 
+  UnDim) 

600  continue 

return 

end 
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subroutine  getlioe(  line) 

Prompting  needs  to  be  done  before  calling  getline. 
Input : 

line  initially  contains  the  default  value 
Return  value: 

line  contains  the  output 

character*80  line*  inline,  blank 
character*!  comment 

data  blank/  '  '/ 
data  comment/ 'CV 

format('  Default:  ',a80) 
f  orma  t (  a  80 ) 


u  itaei 

1020 
Li  endif 


write!  *»  lOOOj  line 
cont inue 

return  nere  if  the  line  is  a  cocmient 
read!  *,  1010,  ena=200)  inline 

DEBUG 

write!  *,  1020)  inline 
format!’  get  line:  '.aSO) 


if!  inline!  1 : 1 )  .eq.cocuaent )  goto  100 
if!  inline. ne. blank)  tnen 
line=inline 

else 

don't  cnange  me  iine....<>de 

endif 

return 

cont inue 

EOF  was  detected 
line  =  blank 


t^deiauit  was  acceptec, 


real  function  getreal!  prompt/ 

character*80  prompt 
real  value 

call  getlinel  prompt) 

read!  prompt,  '(flO.3)'}  value 

„etreai  =  value 
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logical  fuactioa  Agaia(  ierr) 
c 

c  Description:  Print  out  a  diagnostic  and  figure 

c  out  whether  the  user  wishes  to  stop  the  program 

c  or  rerun  the  section  which  just  generated  the  error, 

c 

c  Author:  Chuck  Gendrich 

c 

c  History: 

c  August.  1987.  vl.O 

c 

character*!  YorN 
c 

Li  include  ’’vorticity.e'* 
c 

1000  format(a80/'  Do  you  wish  to  repeat  the  section  where 
*  'the  error  occurred?  [y]') 

1010  format!  alj 

1020  format!'  Please  input  either  "y"'  or  ''n''.') 

c 

10  continue 

write!  *.  1000)  err^  ierr) 
read!  *.  1010.end=20)  YorN 

if!  YorN. ne. 'y' .and. YorU.ne. 'Y* .ana. YorN. ne. '  '.and. 

+  YorN. ne. 'n' .and. YorN. ne. 'N' )  then 
c  the  user  goofed,  can't  type  or  something 

write!  *»  1020) 
goto  10 

else 

Again  =  !  YorN. eq. 'y' .or. YorN. eq. 'Y' .or. YorN. 
return 

end  if 

c 

20  continue 

c 

c  EOF  encountered 

Again  =  .FALSE, 
return 
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subroutine  CutUerev  title) 
c 

Li  include  “vorticity.n"” 
c 

write(  *,  1000)  title 
c 

1000  tori!iat(/ 

0_  _'/ 

x'  ====z  =  =  =  _><_ - cut  here  ror:  '  ,u40/ 

t'  OV) 

c 

return 

ena 
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real  function  d(  VxVy,  dd,  p.  vel,  ierr) 
c 

c  Description:  Take  the  derivative  of  U  or  V  with  respect  to 

c  X  or  y  (as  inaicated  by  dd).  The  derivative  is  to  be  taken 

c  at  point  p  in  velocity  descriptor  array  vel. 

c 

c  du/dy  Procedure:  Look  down.  If  a  'down*  point  exists,  see 

c  what  angle  it's  at.  interpolate  to  find  u  at  a  point  directly 

c  below  the  original  point.  This  tcakes  dx  ==  0. 

c 

c  Then  calculate  du/dy  =  (ul-u2)  /  (yl-y2)  <J  the  point  (xO,  avg(y)) 

c 

c  du/dx  Procedure:  Similar  to  du/dy,  but  swap  X's  and  Y's 

C 

c  dv/dy  Procedure:  Identical  to  cu/dy,  but  use  vel(  Vy)  instead, 

c  This  should  be  taken  care  of  by  the  routine  being  called  with 

c  VxVy  sec  to  either  Vx  or  Vy  appropriate;. 


c  Author:  Chuck  Genaricn 


t  History: 

t  August,  19  87  vl  .0 

c 

ccccccccclccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc 

c 

cLi  aefine  DEBUG 

ct^  aefine  UGllBUG  <— -  only  if  you  want  the  really  gruesome  diagnostics 
a.  deiine  STRICT 

c  only  use  values  corrected  for  ux/dy  (only  use  interpolateu 

c  values) 

Lu  inciuae  ""vorticity .  h” 
c 

integer  VxVy,  ud 
real  oelta 

parameter (  delta  =  39.0; 

real  xO,  x2,  dx 

real  yO,  y2,  uy 

real  angle 


integer  look 

real  interp 

external  look,  interp 

c 

logical  UnDim 

c  --  from  convert...  not  used  nere  - 

^  ifaei  UGHBUG 

a  =  h'umP  t  s  (  vel,  4) 

write(  *,  1000)  1 i, vel( i,X) ,vel( i,Y ) ,vel( i, Vx) ,vel( i, Vy ; , i=l ,n ) 
luOO  iormatC/'  Here  are  cue  velocity  components  and  tr.eir’, 
locations:'/'  X  Y  Vx  Vy'/ 

r  ax,i2,':  ',4(lx,i7.2;;) 

w  end if 
c 

ccccccccclccccccccc2cccccccccaccccccccc4ccccccccc5ccccccccc6ccccccccc7cc 


o  a 


I 
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ierr  =  0 

i£(  dd.eq.ddy)  Chen 

j  =  look(  p,vel.4>  DOWN,  delta,  laaCchs,  ierr; 

Q  i£de£  DEBUG 

write(  *,  1010)  'down',  p,  j 
c*  write(  *,  1010)  'down',  p,  j,  natchs 

LJ  end  if 

else 

j  =  look(  p,vel,4,  RIGHT,  uelca,  'uiatchs,  ierr) 

U  ifdef  DEBUG 

write(  *,  1010)  'right',  p,  j 

*  write!  *,  ,010)  'right',  p,  j,  uatchs 

end  if 
end  if 

1010  rormatC'  "look''d  ',a5,'  from  ',i2,'  ana  found  ’.i2 
c*  T  /'  Here*'s  the  raatch  array; '  ,20Clx,i2j) 
c 

c  error  check  point  ^ 

c 

i£(  J  .eq.BiVD.or .  lerr.ne.O^  then 
d  =  BAD 

if!  lerr.ne.Oy  then 
return 

else 

ierr  =  21 
return 


enaif 

end  if 


>:0 

= 

vel! 

P,a) 

= 

veil 

j  »x) 

ux 

= 

x2  - 

xO 

yO 

vel! 

P ,  i  > 

y2 

= 

vel! 

J 

ay 

yl  - 

>0 

angle  =  canl C  dy,  ax) 
if!  ad.eq.ady)  then 

if!  angle. gt.DOV/K)  then 

k  =  look'.  J  ,  vel  ,4,LLFT,  delta,  :.a tens ,  ierr ; 

else 

k  =  look!  j, vel, 4, RIGHT,  aelta,  i.atchs,ierr 

enaif 

else 

if!  angle. gt. delta)  then 
c  It's  in  tne  4th  quadrant. 

k  =  look!  J, vel, 4, UP,  aelta,  satchs,ierr) 

else 

k  =  1.00X1.  J  ,vfel,4»DOl.’H,  uelta,  ..a tens ,  ierr ; 

enaif 
end  if 


interpolate  li  we  icuna  a  point  to  interpolate  witn. 
if!  it.ne.BAD.ana.  lerr.eq.U)  tnen 
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c  there  is  a  point  to  interpolate  to 

Lj  ifdef  DEBUG 

writeC  *.  1020)  j.k. 

1020  formatC'  Going  to  interpolate  between  points  '.i2i 

+  '  and  '.i2,'.') 

J  end if 

i£(  dd.eq.ddy)  then 

u=interp(  x2,  vel(j,VxVy),  veKk.X).  vel(k.VxVy),  xO.ierr) 
else 

u=interp(  y2,  veKj.VxVy),  veldt, Y),  veldt,  VxVy) ,  y0,ierr) 
endif 
else 

c  use  the  lower/right  velocity  and  ignore  dx/dy...  eeek! 

u  =  vel(  j.VxVy) 
endif 

c 

Lj  ifdet  STRICT 

if(  k.eq.BAD.or. lerr.ne.O;  then 
return 

endif 
o  endif 

c 

du  =  u  -  vel(  p,VxVy) 
if(  dd.eq.ddy;  then 
d  =  du/dy 

call  converC(  xO,  avg(  y0,y2),  du/dy,  0,  UnDim) 
else 

d  =  du/dx 

call  converts  avg(  xO,  x2),  yO,  du/dx,  0,  UnDici) 
endif 
return 

forciat(  3(2x,fS.3}) 


1040 
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subroutine  convert!  xpix,  ypix.  ddpix.  uvpixt  UnDim) 
c 

c  Description:  print  out  the  conversion  to  ’real  units" 

c  of  the  above  values.  Only  print  out  the  ones  which 

c  aren't  zero  (generally  either  ddpix  or  uvpix)... 

c 

c  Output:  all  output  will  be  written  to  stdout.  If  you  want  it  in 

c  a  filet  put  a  tee  on  the  process  when  you  run  it. 

c  e.g.:  lori  38>  vorticity  I  tee  output_file 

c 

c  Procedure:  First  get  conversion  factors  and  the  time  between 

c  each  frame.  Since  this  subroutine  is  ‘‘saved*’,  these  values 

c  only  have  to  be  obtained  once,  then  they're  applied  to  every 

c  subsequent  value  as  appropriate... 


c  Caveats:  Ko  values  are  printec  out  the  first  time  tnrougiu 

c 

c  Author:  ChucA  C-eiicrich 

c 

c  History: 

c  August,  1987  vi .0 


c 

ccccccccclccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc 

c 

cu  aefine  DEBUG 
c 


c 

c 

c 


c 


logical  UnCim 

true  if  values  shoula  oe  non-dimensionalizeu 

real  xpix,  ypix,  udpix,  uvpix 

xpix  ana  ypix  are  (X,Y}  in  pixel  values 
ddpix  IS  a  spatial  aerivative  of  some  velocity 
(e.g.  uu/dx; 

uvpix  IS  tne  product  oi  two  velocities 
(e.g.  Vx*Vy; 

real  xmm,  ^rnsn,  udmm,  uvmm 

X,  y,  dd,  and  uv  converted  to  ‘‘real  units”  or  non- 
dimensiona lized  (local  copies  of  the  numoers  so  that 
we  don't  accidentally  try  to  set  some  constant  equal 
to  something  else..... 

real  HmPixel,  ut.  hmdt,  Xmdt2,  xorf,  yoff 


real  Uw,  nu,  t,  xyract,  adfacc,  uvfact 


c 

Uw 

--  wail  velocity 

- 

nu 

--  kinematic  viscosity 

a 

t 

--  total  elapsed  time 

C 

xy tact 

--  non-dimensionaiizing 

tactor 

for 

K  *  s  unc  y  '  3 

C 

ddf act 

--  non-dimensionalizing 

tactor 

for 

dd '  s 

uvfact 

--  uon-d imensiona 1 izing 

factor 

for 

uv '  s 

real  <^etreal 
external  ^etreai 

cnaracter»60  line 


C 
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aata  limPixelt  at/  2^0.0/ 

format('  Please  enter  the  multiplication  factor  to  convert'/ 
pixels  to  ran.  (pixel_value  *  factor)  =  (imn_value ) ' ) 
fortnat('  Please  enter  the  time  between  frames  (dt).') 

I  format{'  Please  enter  the  Y  offset  for  frame  1.’/ 

+'  (  Y(framel)  -  Yoffset)  =  Y  (abs  dist  to  the  walll’) 
i  format('  V’ould  you  like  the  values  to  be  non-dimensionalized? ' , 

Cn]'} 

I  £armat(a80} 

>  formate '  Please  enter  the  wail  velocity,  Uw.’) 

)  formate '  Please  enter  the  kinematic  viscosity,  nu.') 

I  formate*  Please  enter  the  total  elapsed  time,  t.') 

I  formate/*  convert;  HmPixei;  ',fl0.9.*,  dt :  ',fl0.9,'  Yoff: 

+  flO.7) 

I  rormatt*  convert;  non-dimensional is ing.  nyfact:  ',e9.4, 

adfact:  ',e9.4,'  uviact:  ',e9.4/) 

)  formate*  convert:  printing  values  with  real  units.'/) 

)  formate/*  convert:  Pixel  values  --  l * , f8 .2, ' , ' ,f 8 .2 , ' ) '/ 
spat,  aeriv;  ',fl2.3,'  and  uv  value;  ',fl2.3// 

)  £ormate4e2x,fll.3;) 

;-cccclcccc-cccc2cccc-cccc3cccc-cccc4cccc-cccc5cccc-cccc6cccc-cccc7cc 

get  tne  conversion  to  ram's  and  the  time  between  frames  if  necessary 

If ei^mPixel . nc .0 .and .dt . ne. 0 )  goto  100 

neither  of  these  values  is  permitted  to  be  zero... 

continue 

writel  *,  1000/ 

line  =  '0.029520(309  ;.m/pixei* 

0.0277  ram/pixei  wnen  .50'  =  458.13  pixels 
hmPixei  =  getreaK  Line; 

write!  *,  1020) 

line  =  '0.02000000  sec' 

at  =  getreaiv  line; 

Cradt  =  HmPixei/dt 
;iradt2  =  iimat  *  .Icdt 

write!  *,  1030; 

line  =  ‘560.0  Pixels' 

Yoff  =  getreaiC  Line; 

horf  -  0.0 

cont inue 


,et  non-aimensionalizing  constants  \,or  set  them  to  i.O; 


writei  *,  1040; 
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read(  1050)  line 

if(  lined :  1 )  .eq. 'Y' .or .  lined :  1 )  .eq.  *y’ )  then 
UoDim  =  .TRUE. 

else  i£(  lined :  1 )  .eq. 'U' .or .  lined :  1 )  .eq. ’n’ .or . 

+  lined :  1 ;  .eq. '  ')  then 
UnDim  =  .FALSE. 

else 

goto  20 

end  if 

c 

if(  UnDim)  then 

write!  *.  1060) 
line  =  '5.0 
Uw  =  getreal! 
c 

write!  *,  1070) 
line  =  '0.00310460  in  2/sec.' 
nu  =  ^ecreaiv  iiaey 

C 

write!  *.  1080, 
line  =  ‘5 .25  sec. ' 

t  -  ^etreal!  line; 
c 

xy  tcLcc 
ddfact 
uvfact 
c 

else 

xyfact  =  1.0 
ddfact  =  1.0 
uvfact  =1.0 

end  if 

c 

^  ifoer  DEBUG 

writer  *,  1090;  ;imPi:iei,  at,  Yoif 
if!  UnDim)  then 

write!  *,  1100)  xyfact,  adfact,  uvfact 

else 

write!  *.  1110) 

end  if 
enaif 

c 

if  iL'niPixei  .  eq  .0  .or  .dt  .eq  .0  ;  goto  10 
return 

c  don't  print  anything  out  the  first  time  through 

c 

100  continue 

i_  ifdef  DEBUG 

write;  *,  1120;  xpix,  /pix,  adpix,  uvpix 
_  endif 

xc-.m  =  !xpix  -  Aoff)  *  iimPixei  *  xyfact 
yr-zu  =  !>pi-x  -  Yoff;  *  limPixei  *  xyfact 


=  1.0/!  2.0  *  25.4  *•  sqrtt  nu  *  t)) 

=  4.0  *  sqrt!  nu  *  t)  /  Uw 
=  1.0 

don't  Imow  how  to  non-dimensionalize  this  one 


in./ sec ' 
line) 


aamm  =  ddpix  /  at  ^  adfact 
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c 


c 


uvTom  =  uvpix  *  Mmdt2  *  uvfact 


if( 

ddsnin. 

ne.O.and. 

uvmm.ne.O) 

then 

write( 

*,  1130) 

xmmt  ymm* 

ddmm, 

else 

i£(  ddmm. 

ne.O.and. 

uvmm. eq.O) 

then 

writeC 

*,  1130) 

yirnn,  ymm. 

ddmm 

else 

i£(  ddmm. 

eq.O.and. 

uvmm.ne.O) 

then 

write( 

*.  1130) 

xmm.  ycmit 

uvmm 

else 

write( 

*.  1130) 

xznza.  yiam 

endif 


return 

end 


uvmni 
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subroutine  integrtC  pi.  p2t  vel.  area,  intgrl) 

c 

c  Description:  Integrating  fron  point  pi  to  point  p2  (whose  X-  an& 

c  Y-coordinates  are  described  in  vel).  return  the  area  under  the 

c  curve  and  the  value  of  V  *  ds .  Simple  trigonometric  and 

c  calculus  identities  are  used;  e.g..  the  area  between  a  line  and 

c  the  x-axis  is  1/2  (yl  +  y2)  (x2  -  xi)  and 

c  _  _ 

c  V*ds  =  |V|*ls|*cos(  the  angle  between  the  two  vectors) 

c 

c  Temporary  variables  a/:e  used  to  store  intermediate  results  so 

c  that  the  steps  are  clear, 

c 

c  Author:  Chuck  Gendrich 

c 

c  History: 

c  August,  1987  vl .0 

c 

ccccccccclccccccccc2ccci,ccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc 

c 

cL  define  DEBUG 
u  include  "vor ticity .h" 
c 

integer  pl.  p2 
real  area,  intgrl 
c 

real  avg.  dist.  tanl 
external  avg,  dist,  tanl 
c 

real  Ubar,  Vbar,  V,  s,  thetaV,  thetaS 
c  Ubar  is  the  average  Vx  for  a  side 

c  Vbar  is  the  average  Vy  for  a  side 

c  V  IS  the  magnitude  of  the  velocity  vector 

c  s  is  me  length  or  a  side 

c  thetaV  is  tne  angle  V  makes  with  tne  -^x  axis 

c  thetaS  is  the  angle  S  makes  with  the  -^x  axis 

c 

real  pilSO 

parameteri  piloO  =  0.017453293) 
c  pi/ 180  for  converting  degrees  to  radians 

c 

if(  pl  .eq.O  .or . p2  .eq.O)  then 

c  we're  integrating  to  or  from  a  non-existent  point 

area  =  0.0 

intgrl  -  0.0 

return 

end  if 

xl  =  vel(  pl,  X) 

yi  -  vel(  pl,  Y) 

c 

x2  =  vei(  p2, 

y2  =  vel(  p2,  Y; 

c 

GX  =  x2  -  xi 

ay  -  y2  -  yi 
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c 

Ubar  =  vel(  pi,  Vx) 

Vbar  =  vel(  pi,  Vy) 

c  Ubar  =  avg(  vel(pl,Vx),  vel(p2,Vx)) 

c  Vbar  =  avg(  veUpl.Vy),  vei(p2,Vy)) 

c 

V  =  sqrt(  Ubar*Ubar  ~  Vbar*Vbar) 
s  =  sqrt(  QX*dx  +  dy*dy) 
c 

thetaV  =  canK  Vbar,  Ubar) 
thetaS  =  tanl(  dy,  dx) 
c 

intgrl  =  V  *  s  *  cos(  (thetaV  -  tfaetaS)  *  pil80) 
area  =  0.5  *  (yl  ^  y2)  *  dx 
c 

U  ifdef  DEBUG 

write(  *•  1000)  pi,  p2,  dx,  uy,  Ubar,  Vbar,  V,  s,  thetaV,  chetaS 
1000  formate/’  intgrt;  from  ',i2,'  to  ',i2,-.  ax:  ',£8.3,'  dy: 

T  £8.3/’  Ubar;  ',£8.3,'  Vbar;  '.£8.3,'  ->  V:  ',£8.3/ 

's:  ',£8.3,'  thetaV:  ',£9.3,'  and  thetaS;  ',£9.3) 
write(*,  1010)  intgrl.  area 

1010  foriaatC'  line  intgri:  ',£12.3,'  urea:  ',el4.3/) 

U  end if 
c 

return 

end 


I 
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subroutine  vector(  dy»  dx.  scale) 

c 

c  Description:  Draw  a  2D  vector  with  arrow  head  from  the  current 

c  point  through  the  final  point.  The  length  of  the  vector  will  be 

c  equal  to  its  original  length  times  "scale.** 

c 

c  The  arrow  head  will  always  be  the  same  size  (more  or  less), 

c 

c  Current  values  for  linewidth.  color,  etc.  are  used, 

c 

c  Author:  Chuck  Gendrich  <cpg> 

c  ("*You  mean  I  can't  just  find  this  in  some  library?  !CC  ) 

c 

c  History: 

c  12  dec  87  —  original  inception  and  coding 

c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 


c 


c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 


real  uy.  dx.  scale 

real  SideLen.  dAng 
parameter (  SideLen  =  8.0; 
parameter (  dAng  =  20.0) 

the  length  of  each  side  of  the  arrow  (in  pixels) 
the  angle  that  each  side  makes  with  the  vector 

real  xl.  yl,  x2.  y2.  x3.  y3 .  a2.  a3 

destination  points  for  the  far  end  of  the  vector  and  the 
two  ends  of  the  arrow  head... 

real  angle 

the  angle  chat  the  vector  makes  with  the  X-axis 

real  can! 
external  tanl 
real  pi 

parametar(  pi  =  3.1A1592653589) 

Figure  the  ant,ies  that  the  lines  will  be  at  VJRT  the  X-axis 

angle  =  tanl ^  uy ,  ax) 

a2  =  (angle  -  uAng)  *  pi  /  180.0 

a3  =  (angle  +  dAng)  *  pi  /  180.0 

Figure  the  relative  moves 


xl  =  scale  *  dx 
yl  =  scale  *  dy 


x2  = 

-SideLen 

* 

cos  ( 

a2) 

y2  = 

-SideLen 

* 

sin( 

a2) 

x3  = 

-SideLen 

★ 

cos  t 

a3) 

y3  = 

-SideLen 

* 

sint 

a3; 

Draw  It 


call  rardt  xl,  yl ) 
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logical  function  GoOnO 

character*!  char 
logical  first,  slow 
save 

data  first/ .TRUE./ 

if(  first)  then 

first  =  .FALSE, 
print*. 'Slow  mode?  [y]' 
read(  *.'(Al)')  char 
if(  char .eq. 'n' .or .char .eq. 'N* 
slow  =  .FALSE. 

else 

slow  =  .TRUE. 

end  if 

end  if 

if(  SLOW)  then 

print*. ' <RETURN>  to  continue.. 
read(  *,'(A1)')  char 

endif 

GoOn  =  .TRUE. 

return 

end 


)  then 
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I 

i 


G 

D  makefile  for  the  vorticity  determination  program  (vorticity) 

U  and  other  miscellaneous  test  and  error-analysis  programs 

0 

LIBS= 

DESTDIR=/usr/thor/gendrich/bin/ iris/ 

0BJS1=  getline.o  getframe.o  errors. o  sort.o  inatch_points .o  count_points .o  \ 
velocities. o  make_poly.o  geometry. o  look.o  main.o  fluids. o  \ 
cut_here.o  ddx.o  convert. o  integrate. o  vector. o  GoOn.o 
0BJS2=  TestLook.o  getline.o  getframe.o  count_points.o  geometry. o  sort.o  \ 
convert. o  GoOn.o 

0BJS3=  TestMatch.o  getline.o  getframe.o  errors. o  look.o  count_points .o  \ 
geometry. o  sort.o  convert. o  vector. o  GoOn.o 
0BJS4=  TestDdx.o  getline.o  getframe.o  errors. o  sort.o  inatch_points .o  \ 
count_points.o  velocities. o  geometry. o  look.o  ddx.o  convert. o  \ 
vector. o  GoOn.o 

0BJS5=  getline.o  getframe.o  errors. o  sort.o  match_points .o  count_points.o  \ 
velocities. o  one_poly.o  geometry. o  look.o  TestPoly.o  fluids. o  \ 
cut_here.o  ddx.o  convert. o  integrate. o  vector. o  GoOn.o 
0BJS6=  Errors. o  stats. o  vector. o  sort.o  count_points .o  match_points .o  \ 
getframe.o  getline.o  geometry. o  look.o  errors. o  convert. o 
OBJS=  S(OBJSl)  $(0BJS2)  S(0BJS3)  $(0BJS4)  one_poly.o  $(0BJS6) 

SRCS1=  main.f  getframe.f  count_points.f  sort.f  match_points.f  velocities.!  \ 
geometry.!  make_poly.f  look.!  fluids.!  getline.f  errors.!  cut_here.f  \ 
ddx.f  convert. f  integrate.!  vector. f  GoOn.f 
SRCS2=  TestLook.f  getline.f  getframe.f  count_points . f  geometry. f  sort.f 
SRCS3=  TestMatch.f  look.f  getframe.f  errors.!  getline.f  count_points . f 
SRCS4=  TestDdx.f  ddx.f  getframe.f  count_points . f  sort.f  match_points . f  \ 
velocities.!  geometry.!  look.f  convert.!  getline.f  errors. f 
SRCS6=  Errors.!  stats.! 

ALL=  vorticity  TestLook  TestMatch  TestDdx  TestPoly  Errors 

0= 

LI  can't  optimize  FORTRAN  cone,  bummer. 

DBG= 

LiDBG=  -g  -DDEBUG 
C0PTS=  -Zf  -Zg 

CFLAGS=  oO  S(COPTS)  $(CC0KFIG)  >i(DBG) 

CC=  f77 

INSTALL=n;v  -f 
CHM0D=  chmoc 
OBJMODE=755 

.  f  .0 : 

$(CC)  SCCFLAGS)  -c  -o  S*.o  $*.f 
vorticity:  $(0BJS1) 

S(CC)  l?(COPTS)  -o  vorticity  $(0BJS1}  3(LIBS) 

$(CHI'j'OD)  $(0BJ1X)DE)  vorticity 

TestLook:  S(0BJS2/ 

;>(CC;  $(COPTS)  -o  TestLook  5>(OBJS2}  ;?CLIBS) 
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$(CHH0D)  $(OBJMODE)  TestLook 

TestMatch:  $(0BJS3) 

$(CC)  $(COPTS)  -o  TestMatch  $(0BJS3)  $(LIBS) 
$(C11M0D)  $(0BJM0DE}  TestMatch 

TestDdx:  $(0BJS4} 

$(CC)  $(COPTS)  -o  TestDdx  $(0BJS4)  $(LIBS) 
$(CHMOD)  $(OBJMODE)  TestDdx 

TestPoly:  $(0BJS5) 

$(CC)  $(COPTS)  -o  TestPoly  $(0BJS5)  $(LIBS) 
$(CHMOD)  $(0BJM0DE)  TestPoly 

Errors:  $(0BJS6) 

$(CC)  $(COPTS)  -o  Errors  $(0BJS6)  $(LIBS) 
$(CHMOD)  $(0BJM0DE)  Errors 


install;  vorticity  TestPoly  Errors 

$(INSTALL)  vorticity  S(DESTDIR) 
$(INSTALL)  TestPoly  $(DESTDIR) 
S(INSTALL)  Errors  $(DESTD1R) 

all:  p(ALLJ 


print : 

pr  $(SRCS1)  Makefile  I  rsh  freyja  Ip 


clean ; 

nn  -f  core  junk  S(ALL)  $(OBJS) 
mv  vor<*  Results 


TestLook. 0 : 

TestLook. f  vorticity. h  graphics. h 

Testtuitch.o : 

TestMatch. f  vorticity. h  graphics. n 

TestDdx. o: 

TestDdx. £  vorticity. h  graphics. h 

TestPoly .0 : 

TestPoly. f  vorticity. h  graphics. h  Makefile 

main.o : 

main.f  vorticity. h  graphics. h  Makefile 

Errors .0 : 

Errors.!  stats. h 

stats. o: 

scats. f  stats. n 

convert .o : 

c  onver  t . £  Make! i 1 e 

count_points .o: 

count_points . £  vorticity. n 

cut_here .0 : 

cut_here.£  vorticity. h 

ddx.o: 

ddx.f  vorticity. h  Makefile 

errors .0 : 

errors. f  vorticity. e 

f luias.o: 

fluids.!  vorticity. h 

geometry .0 ; 

geometry.!  vorticity. h  Makefile 

getline.o: 

get  line.!  Makefile 

getf  raise,  o: 

getframe.f  vorticity. h 

integrate .0 : 

integrate.!  vorticity. h  Makefile 

look.u ; 

look.!  vorticity. h 

make_poly .0 : 

make_poly.f  vorticity. h 

one_poly .o : 

one_poly.!  vorticity. h 

.aatch_points  .0 : 

inatch_points . !  vorticity. h 

sort .0 : 

sort.!  vorticity. h 

velocities . o : 

velocities.!  vorticity. h 
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ABSTRACr 

A  non-intrusive  photochemical  phenomena  called 
photochromism  is  used  to  obtain  important  fluid  quantities 
such  as  velocity,  and  vorticity  over  a  two-dimensional  domain 
of  a  fluid  flow.  Using  a  uv  laser,  and  a  novel  beam  splitting 
and  steering  technique,  a  grid  of  lines  is  impressed  in  the 
prepared  fluid.  Two  successive  pictures  are  needed  to  calculate 
the  fluid  dynamic  quantities.  Measurements  in  a  Stokes'  layer 
are  compared  against  exact  solutions,  and  measurement  of  the 
circulation  of  a  vortex  ring  is  presented.  The  technique  will  be 
useful  in  turbulent  and  unsteady  flows. 

1.  INTRODUCTION 

There  is  an  increasing  need  to  obtain  flow  field 
information  over  a  region  of  flow  simultaneously.  Examples 
include  the  study  of  mixing  in  internal  combustion  engines,  and 
the  understanding  of  the  importance  of  coherent  motions  in 
many  turbulent  and  unsteady  flows.  Furthermore,  there  is  a 
need  to  obtain  more  detailed  information  about  these  and 
many  simpler  flows,  including  vorticity  and  strain-rate 
measurements.  Existing  probe  techniques  are  at  best  able  to 
provide  single  point  estimates  of  these  quantities  in  flows  with 
steady  mean  velocities. 

Hummel  and  his^ co-workers  first  exploited  the  phenomena 
of  photochromism  ’  .  Using  a  single  line  of  uv  radiation  for 
excitation,  he  and  his  co-workers  were  able  to  measure  one 
component  of  velociw  in  flows  with  predominately  one 
direction  of  motion.  With  the  advent  of  high  power  uv  lasers 
and  new  chemicals  it  is  increasingly  possible  to  apply  the 
technique  to  more  complicated  flows.  We  have  generalized  the 
technique  to  enable  it  to  uniquely  tag  specific  points  in  the  flow 
on  a  specified  grid.  As  a  result  we  have  a  non-intrusive  photo- 
optical  technique  which  can  be  used  to  obtain  two  velocity 
components  and  the  spatial  gradients  of  these  components,  at  a 
large  number  of  points  in  a  flow  field. 

2.  USE  OF  PHOTOCHROMISM  TO  MEASURE  FLUID 
DYNAMIC  QUANTITIES 

Thp  phenomena  of  photochromism  has  been  defined  by 
Brown’  as  "a  reversible  change  of  a  single  chemical  species 
between  two  states  having  distinguishably  different  ab.sorption 
spectra,  such  change  being  induced  in  at  least  one  direction  by 
the  action  of  electromagnetic  radiation."  The  chemical  used  in 
this  study  was  Kodak  l,3-Trimethyl-8-nitrospiro[2-H-l- 
benzopyran-2,2’-indoline].  We  excited  it  at  351nm  using  pulses 
of  20ns  produced  by  a  Tachisto  Excimer  laser  running  on  XeF 
(approximately  lOOmj  output).  Dissolved  in  the  working  fluid, 
which  was  kerosene,  we  found  that  the  chemical,  thus  excited, 
reflected  blue  light  for  the  order  of  tens  of  seconds.  Thus,  the 
fluid  was  tagged  long  enough  for  it  to  move  appreciably.  In 
practice,  for  the  velocities  of  our  experiments,  we  needed  to 
follow  the  lines  only  a  fraction  of  a  second  to  obtain  the 
information  required. 


2.1  Use  of  a  grid  to  tag  fluid  particles 

To  enable  us  to  measure  two  components  of  the  velocity  in 
a  flow  using  the  visualization,  we  must  tag  and  follow  specific 
fluid  particles.  The  simplest  way  to  do  this  is  have  two  lines 
intersect.  That  intersection  defines  a  point  in  the  flow  which, 
because  it  is  colored,  we  can  follow.  To  measure  a  spatial 

fradient,  we  must  simultaneously  mark  two  points  in  the  flow. 

o  enable  us  to  measure  a  spatial  gradient  in  a  prescribed 
direction,  we  need  several  points  and  must  use  interpolation. 
Thus,  creation  of  a  grid  of  intersecting  marked  fluid  lines 
would  enable  both  velocities  and  gradients  to  be  measured. 


2,2  The  beam  divider 

We  sought  a  method  of  dividing  a  beam  of  laser  light  into 
’n’  beams  where  ’n’  would  be  5-25,  using  a  static  device.  In 
doing  so,  we  sought  to  lose  as  little  light  as  possible.  Thus,  we 
designed  a  divider  that  appears  as  an  oversized  blazed 
reflection  grating;  the  steps  size  ’a’  is  approximately  2mm.  The 
incident  laser  ener©i  was  directed  at  a  shallow  angle  to  the 
plane  of  the  grating.  Becau.se  the  facets  are  so  much  larger 
than  a  wavelength  of  light,  we  expected  to  primarily  see 
specular  reflection.  The  separation  of  the  specularly  reflected 
portions  of  the  beam  were  determined  by  the  difference  of  the 
angle  of  incidence  and  the  blaze  angle,  and  the  grating  step 
width  and  spacing  as  shown  in  Fig.  1.  A  photo  of  the  actual 
beam  divider  is  shown  in  Fig.  2.  It  was  fabricated  by 
individually  fixing  Aluminum  coated  mirrors  (with  the  coating 
optimized  for  shallow  incidence  angles)  to  a  machined  base. 
Apart  from  diffraction  effects,  no  incidence  laser  energy  would 
be  lost. 


Fig.  1.  The  specular  reflection  from  the  beam  divider. 


As  expected,  however,  diffraction  effects  were  exhibited. 
A  photo  of  the  output  of  the  beam  divider  using  a  Helium 
Neon  laser  is  shown  in  Fig.  3.  We  .see  the  specularly  reflected 
beams,  and  the  superimposed  diffraction  pattern.  The 
diffraction  bands  are  spaced  as  expected  from  the  grating 
equation*.  The  appearance  of  the  first  and  higher  order  bands 
are  a  source  of  ’noise’  in  the  present  grid  making  .scheme. 
When  experiments  were  performed  with  the  Excimer  laser  at 
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Fig.  3.  A  photo  of  the  output  of  the  beam  divider. 

Fig.  4  shows  a  sketch  of  the  optical  set-up  for  making  the 
grid  shown  in  Fig.  5a.  Note  that  the  lines  in  Fig.  5a  remain 
uniform  and  show  little  divergence.  The  mesh  size  is  the  order 
of  Itnm.  The  two  photos  in  Fig.  5  are  0.4  seconds  apart,  taken 
by  a  still  camera  with  a  motor  drive.  From  the  distortion  of  the 
grid  in  Fig.  5b  with  respect  to  that  in  Fig.  5a,  (one  laser  pulse, 
two  photos  At  apart)  we  can  calculate  all  the  kinematic 
quantities  of  the  flow  in  the  plane  of  the  film.  (Note,  that  in 
practice  we  would  not  use  a  second  image  which  has  distorted 
as  much  as  Fig.  5b,  but  we  present  it  to  illustrate  the  ideas 
involved.) 


(b) 


Fig.  5.  An  example  of  the  grid,  a)  shortly  after  the  laser  pulse, 
b)  .4sec  later. 

2.3  Algorithms  to  obtain  fluid  dynamic  quantities 

For  convenience,  we  pick  up  one  ’grid  box’  in  Fig.  5a 
which,  if  the  mesh  size  is  small  enough,  can  be  thought  of  as  a 
text  book  fluid  ’particle’  in  flow  at  time  t  =  0.  ^is  fluid 
particle  is  shown  using  solid  lines  in  Fig.  6a.  As  this  particle 
moves  with  the  flow,  it  may  undergo  several  motions  such  as 
translation,  rotation  and  deformation.  The  dashed  lines  in  Fig. 
6b  show  the  same  fluid  particle  after  a  shoit  time  interval.  At. 
Because  we  know  the  history  of  this  specific  fluid  particle,  we 
can  compute  the  velocity  of  each  comer  by  taking  the  time 
derivative  of  displacement  ds  associated  with  each  point,  that 
is,  ds/At.  The  deduction  procedures  are  shown  in  Fig.  6b. 


SPIE  Vol.  814  Photomechanics  and  Speckle  Metrology  (1987)  /  707 


Note  that  V  is  desimated  as  the  velocity  at  point  1  which  is  To  obtain  velocity  gradients,  vorticity,  and  strain  rates 
then  positioned  at  the  midpoint  between  point  1’  and  point  1".  from  the  velocities,  we  must  difference  across  the  distance  of  a 
Fig.  7  shows  the  velocity  vectors  deduced  from  Fig.  6.  Actually,  grid  box.  Thus,  for  a  given  flow  the  grid  must  be  small  enough 
a  fluid  particle  moving  in  a  general  three-dimensional  flow  to  resolve  the  smallest  scales  of  motion, 
field  may  have  motions  about  all  three  coordinate  axes, 

however,  we  limit  the  discussion  to  the  projection  of  this  An  alternative  method  can  be  used  if  we  are  only 

motion  onto  the  plane  of  the  photos  which  is  parallel  to  the  interested  in  the  vorticity.  Using  the  definition  of  the 
initial  plane  of  the  grid.  circulation  r. 
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r  -  /V-ds  (1) 
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and  Gauss’  theorem  to  relate  the  surface  integral  to  an  area 
integral. 


r  =  f  w-ndA 

A 

we  can  obtain  the  vorticity  component  normal  to  the  film 
plane,  from  a  grid  box  of  arbitrary  shape,  with 
circumference  C  and  area  A. 

The  velocity  component  V-ds  in  Eq.  (1)  is  estimated  by 
forming  the  average  of  the  corner  velocity  components  (see 
Fig.  7).  Dividing  r  by  the  area  of  the  fluid  particle  results  in 
the  average  vorticity  at  the  centroid  of  this  fluid  particle.  Tliis 
technique  has  the  advantages  of  avoiding  a  second  differencing 
of  the  experimental  data. 

Thus,  by  following  the  distortions  of  a  single  grid  box  we 
can  obtain  the  velocities,  circulation,  vorticity  and  even  the 
Reynolds  stresses  and  strain  rates  if  desired,  over  that  small 
area.  By  doing  this  for  the  other  grid  mesh  points,  we  can 
obtain  the  important  fluid  dynamic  ouantities  at  many 
locations  over  a  two-dimensional  field  in  a  fluid  flow 
simultaneously. 

2.4  Error  analysis 

Consider  the  grid  box  in  Fig.  6.  We  will  assume  that  the 
grid  mesh  is  spaced  <  and  that  a  worse  case  line  movement 
betweei^  images  is  1/2  the  mesh  width.  After  time  At  standard 
analysis  shows  that  the  error  in  estimating  the  vorticity  is 


Fig.  6.  Procedure  for  conversion  of  the  displacement  of  a  grid 
box  to  velocities. 


to  first  order  in  e. 

As  an  example  of  the  designed  accuracy  of  our 
experiments,  if  we  produce  lOO/im  lines  on  a  grid  mesh  of  Itnm, 
and  our  error  in  reading  the  center  of  the  Tines  is  10%,  we 
obtain  ««/«  =  .01.  For  At  =  .Olsec,  the  error  in  the 

measurement  of  vorticity  is  1/sec.  This  is  better  than  has  been 
possible  with  hot-wire  probes,  and  is  comparable  with  (jhe  most 
sophisticated  single  point  laser  doppler  measurements. 

3.  EXPERIMENTAL  RESULTS 

Two  experiments  were  performed  to  demonstrate  the 
technictue.  Both  used  deodorized  kerosene  with  10  ppm  of  the 
photocnromic  chemical  dissolved  in  it.  The  data  were  recorded 
on  Kodak  Tmax  400  ASA  film  through  a  Photo-Sonics  35mm 
high  speed  movie  camera,  framing  between  150  and  200  frames 
per  second.  The  films  were  analyzed  on  a  NEC  film  analyzer 
which  is  accurate  to  .05mm. 


Fig.  7.  Decomposition  of  the  velocities  used  to  calculate  the 
circulation  around  a  grid  box. 
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We  can  test  the  accuracy  of  the  techniaue  in  a  flow  for 
which  an  exact  solution  of  the  Navier-Stokes  equations  is 
obtainable.  The  flow  is  that  of  an  impulsively  started  infinite 
plate  in  an  initially  still  infinite  ambient.  The  problem  involves 
the  two-dimensional  diffusion  of  vorticity  generated  at  the 
plate  into  the  ambient.  The  thickness  of  the  boundary  layer,  s, 
IS  a  function  of  the  kinematic  viscosi^,  u,  and  the  time  after  the 
start  of  the  plate’s  motion.  We  have  approximated  this  motion 
using  a  belt  in  a  tank  which  can  be  rapidly  started.  Fig.  8  shows 
the  experimental  setup.  The  width  of  the  belt  is  17.8  cm,  and 
the  distance  between  the  rollers  is  152.4  cm,  with  the  test 
position  90  cm  downstream  of  the  lead  roller.  All  experimental 
data  was  obtained  before  the  leading  edge  effects  reached  the 
test  position.  Furthermore,  the  aspect  ratio  is  sufficient  to 
prevent  the  disturbances  generated  in  the  corners  from 
reaching  the  center  of  the  belt  at  the  test  section  in  the  test 
time  interval.  Therefore,  the  wall  layer  flow  measured  was  two- 
dimensional.  The  speed  of  the  belt  was  adjustable  from  2.5  to 
22.9  cm/sec.  It  reaches  a  constant  speed  within  a  second.  This 
short  acceleration  reasonably  approximates  an  impulsive  start. 
A  grid  mesh  size  of  .15*  was  used,  with  the  laser  lines 
intersecting  as  in  Fig.  8. 


Fig.  9.  Experimental  measurements  of  the  velocity  gradient  in 
a  Stokes’  layer  compared  wnth  the  theoretical  solution.  A  check 
of  the  continuity  equation  3v/ay  =  0  is  also  shown. 

TTie  experimental  setup  used  to  generate  the  vortex  ring  is 
shown  in  Fig.  10.  A  small  amount  of  fluid  from  a  constant  head 
reservoir  was  emitted  from  a  1.9  cm  orifice  to  create  the  rings. 
A  solenoid  valve  whose  opening  time  could  be  adjusted 
controlled  the  outflow.  The  kerosene  in  the  reservoir  was  dyed 
with  chlorophyll  to  partially  mark  the  ring.  The  convection 
velocity  of  the  rings  was  measured  to  be  3.8  cm/sec  at  the  test 
section  which  is  located  7.6  cm  downstream  of  the  orifice.  The 
Reynolds  number  of  the  ring  is  approximately  700,  based  on 
the  ring’s  core  diameter  and  its  convection  velocity,  ’^e  flow  is 
axisymmetric,  so  that  we  only  needed  to  align  our  grid  through 
an  axis  of  symmetry  and  measure  the  flow  in  half  of  the  ring  to 
measure  its  circulation. 


Results  of  the  velocity  gradient  obtained  from  double 
differentiation  in  the  Stokes  layer  are  shown  in  Fig.  9,  where 
they  are  compared  writh  the  exact  solution.  Everything  is  non- 
dimensionalized  bv  the  similarity  variables  of  the  exact 
solution.  Because  the  flow  is  two-dimensional,  incompressible, 
and  doesn’t  change  with  the  streamwise  coordinate,  x,  the 
continuity  equation  reduces  to  3v/dy  =  0,  and  a  check  on  this  is 
also  shown  in  Fig.  9  (v  is  the  velocity  component  normal  to  the 
wall,  and  y  is  the  coordinate  normal  to  the  wall).  We  see  that 
both  the  error  in  the  estimates  of  the  continuity  equation,  and 
in  the  velocity  gradient  distribution  are  of  similar  magnitude. 
The  absolute  value  of  the  error  is  ±2/sec,  consistent  with  the 
design  error  analysis.  The  deviation  of  the  experimental  profile 
in  the  near  wall  region  is  a  consequence  of  the  grid  mesh  being 
too  large  to  adequately  resolve  the  large  value  of  the  gradient 
au/ay.  On  the  whole  the  results  are  very  encouraging. 
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3.2  Vortex  ring  flow 

The  real  value  of  the  technique  is  its  ability  to  measure 
flows  that  are  not  measurable,  analytically  solvable  or 
numerically  calculable.  The  vortex  ring  is  a  common  place 
example.  It  represents  a  class  of  flow  fields  that  are  of 
particular  importance  in  turbulent  and  unsteady  flows. 


Fig.  10.  Setup  used  to  create  and  measure  the  circulation  of  a 
vortex  ring. 

The  instantaneous  circulation  distribution  of  the  ring  is 
shown  in  Fig.  11.  Although  no  solution  exists  for  the^ 
experimentally  created  ring,  we  show  a  comparison  with  Hill’s 
exact  .solution  for  a  ring  in  which  the  vorticity  occupies  a 
sphere.  The  core  of  our  ring  is  smaller  than  that  of  the 
hypothetical  Hill’s  ring,  the  overall  siiape  is  an  oblate  spheroid, 
and  it  is  unlikely  that  the  vorticity  distribution  is  similar,  but,  as 
the  data  shows,  the  circulation  of  the  ring  (Reynolds  number  = 
700)  is  not  very  different. 
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Fig.  1 1.  The  circulation  distribution  of  a  laminar  vortex  ring 
(Reynolds  number  =  700),  compared  with  Hill’s  solution  for  a 
spherical  vortex. 


4.  DISCUSSION 

The  demonstrated  accuracy  of  the  technique  is  more  than 
adequate  for  many  fluid  dynamic  problems.  Since,  in  principle, 
the  calibration  depends  oivly  on  a  time  and  a  length  scale,  there 
is  further  room  for  improvement.  In  practice,  determining  the 
intersection  points  of  the  grid  lines  is  limited  by  the  grain  size 
of  the  film  used,  and  the  distortion  of  the  developing  process. 
Improvements  in  the  grain  limitation  can  be  made  by  using  a 
large  format  still  camera  with  a  shuttering  mechanism  that  will 
enable  a  double  exposure  of  the  grid  to  be  made  on  a  plate. 
Furthermore,  image  processing  techniques  can  be  used  to 
increase  the  accuracy  and  repeatability  of  the  results.  We  are 
currently  pursuing  these  avenues. 

The  technique  is  not  restricted  to  slowly  varying  flows  as 
are  so  many  flow  visualization  techniques.  Furthermore  it  is 
not  restricted  to  Newtonian  fluids.  We  feel  that  the  whole  field 
approach  can  supplement  LUA  measurements  where  velocity 
fields  are  required,  as  well  as  going  beyond  to  make  Reynolds 
stress,  vorticity,  and  strain  rate  measurements  where  needed. 

Since  the  photochromic  grid  following  technique  requires 
only  a  single  pulse  of  uv  light  to  create  the  grid  in  a  suitably 
mixed  solution,  the  technique  can  be  used  in  many  labs  which 
possess  single  pulse  or  low  pulse  rate  lasers. 

Current  limitations  include  to  need  to  work  in  organic 
liquids,  required  to  dissolve  the  chemicals,  and  available 
penetration  of  the  laser  beam  into  the  fluid.  This  problem, 
addressed  in  Ref.  2  is  a  function  of  both  the  laser  power  and 
the  chemical  concentration.  It  can  be  alleviated  somewhat  by 
scanning  the  laser  beam  across  the  ’steps’  in  the  beam  divider, 
thereby  putting  the  full  energy  of  the  laser  in  every  line.  Since 
the  excitation  takes  place  in  picoseconds,  during  a  20 
nanosecond  pulse  width,  the  scanning  (through  minutes  of  a 
degree)  can  fee  accomplished  with  off  the  self  equipment,  and 
each  line  in  the  fluid  will  be  proportionally  more  intense,  or 
the  chemical  concentration  reduced  to  enable  la.ser 
penetration  further  into  the  fluid. 


5.  CONCLUSIONS 

A  new  technique  for  obtaining  velocity,  and  velocity 
gradient  data  over  a  two-dimensional  domain  of  a  fluid  flow  is 
presented.  Calibration  depends  only  upon  a  time  and  a  length 
scale.  Coinparison  of  measurements  with  the  exact  solution  in 
a  Stokes’  flow  indicates  that  its  accuracy  can  be  predicted  by 
classical  analysis,  and  is  comparable  or  better  than  that 
achievable  with  the  best  single  point  probe  techniques.  The 
circulation  of  a  vortex  ring  has  been  measured. 
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ABSTRACT 


By  dissolving  a  photochromic  chemical  in  an  organic 
liquid  and  exciting  it  with  ultraviolet  light,  a  colored  line 
can  be  created  anywhere  in  the  liauid.  A  device  has 
been  constructed  which  divides  a  laser  beam  into  ’n’ 
lines  which  are  made  to  intersect  to  form  a  md. 
Introducing  two  sets  of  these  lines  into  the  flowing  fluid 
produces  a  grid  of  lOOtim  colored  fluid  lines  with  mesh 
spacing  of  approximately  1mm.  The  color  change 
persists  long  enough  so  that  the  grid  lines  are  distorted 
by  the  flow  before  disappearing.  The  process  is 
completely  non-intrusive.  Tne  lines  are  atudyzed  by  an 
image  processing  algorithm  to  allow  the  grid  intersection 
points  to  be  located  quickly  and  with  high  precision. 
The  movement  of  the  intersection  points  in  a  unit  of 
time  results  in  temporal  gradients  ot  position,  which  in 
conjunction  with  estimates  of  the  circulation  around 
each  grid  box  can  be  used  to  obtain  the  vorticity. 
Additionally,  differencing  again  across  a  grid  box  results 
in  spatial  gradients  of  velocity.  Experiments  in  a  Stokes’ 
flow  indicate  that  the  vorticity  can  be  measured  to 
1/sec.  Application  to  the  measurement  of  the  formation 
of  lift  on  an  impulsively  staned  airfoil  has  been  made 
and  compared  with  known  results.  It  appears  that  the 
separating  streamline  is  fixed  to  the  trailing  edge 
significantly  before  the  Kutta  condition  is  established, 
resulting  in  a  starting  vortex  with  less  circulation  than 
expected. 


1.  INTRODUCTION 


The  generation  of  lift  has  important  technological 
consequences  for  both  military  and  civilian  aircraft. 
Whenever  a  change  in  velocity,  attitude,  or  wing 
conformation  is  made,  lift  changes.  For  helicopters, 
each  blade  on  the  rotor  continually  experiences  velocity 
changes.  Batchelor^  states:  "It  is  a  remarkable  fact  that 
in  praaice  a  circulation  is  generated  around  an  airfoil, 
owing  to  the  convection  of  a  non-zero  amount  of 
vonicity  from  the  rear  edge  of  the  airfoil  at  an  initial 
stage  of  motion,  and  that  when  the  airfoil  is  in  steady 
motion  the  circulation  is  established  with  just  this  specie 
value."  The  subject  of  this  paper  is  to  present 
measurements  using  a  new  experimental  tool  that  can 
provide  the  detailed  data  needed  to  critique  our 
understanding  of  this  problem.  It  enables  the 
simultaneous  measurement  of  fluid  velocity  and  vorticity 
at  a  large  number  of  points  in  a  flow  field.  We 
furthermore  present  preliminary  measurements  that 
suggest  this  problem  is  not  as  well  understood  as  it 
appears,  even  in  the  simplest  example  of  an  impulsiveiv 
started  two-dimensional  wing  (  called  an  airfoil). 

The  predominate  view  of  the  mechanism  for  the 
generation  of  lift  on  an  airfoil  with  a  sharp  trailing  edge, 
is  that  lift  occurs  because  the  zero-lift  irrotational  flow 
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field  which  is  set  up  immediately  after  an  impulsive  start 
would  lead  to  the  physically  unrealistic  condition  of  fluid 
passing  around  the  trailing  edge  which  has  zero 
curvature  at  infinite  velocity.  If  considered  allowable  at 
the  initial  instant,  when  the  viscous  layer  around  the 
airfoil  is  infinitely  thin,  this  fluid  would  move  into  a 
strong  adverse  pressure  gradient  in  the  upper  surffice  of 
the  airfoil  and  separate  because  the  retarded  fluid  in  the 
boundary  layer  on  the  lower  surface  does  not  have  the 
energy  to  negotiate  the  pressure  rise  from  the  infinite 
suction  at  the  trailing  edge  to  the  rear  stagnation  point 
(see  Batchelor^  and  Lighthill^  ^).  A  vortex  forms  and  is 
shed  from  the  body.  When  this  happens,  the  rear 
stagnauon  streamline  moves  to  the  trailing  edge.  Shortly 
after,  a  smooth  flow  is  established,  and  we  have  the 
Kutta  condition,  which  states  that  for  the  conditions  of 
this  problem,  the  separation  streamline  leaves  from  the 
trailing  edge,  and  there  is  no  net  shear  across  this 
streamline.  Figure  1  summarizes  this  classical  view. 
The  shed  vortex  is  the  so  called  ’starting  vortex',  which 
many  authors  suggest  contains  the  full  circulation. 
Calculations  of  the  roll  up  of  a  semi-infinite  vortex  sheet 
such  as  those  of  PuUin"  support  this  point  of  view. 


Figure  1.  a)  Evolution  of  the  starting  vortex  from  the 
irrotational  solution,  b)  The  starting  vortex  after  the 
Kutta  condition  has  been  satisfied. 


Many  authors  have  not  been  entirely  comfonable 
with  this  view.  Lamb^,  states  "It  is  still  not  altogether 
easy,  in  spite  of  the  attempts  that  have  been  made,  to 
trace  out  deductively  the  stages  by  which  the  result  is 
established  when  the  relative  flow  is  started." 
Goldstein  .  after  presenting  it,  states:  "In  this  way  the 
establishment  of  such  a  circulation  may  be  held  to  be  at 
least  partly  explained".  Prandtl^^  states  that  "when  the 
motion  begins,  a  surface  of  discontinuity  forms  at  the 
rear  edge  of  the  airfoil,  and  the  fluid  originally  flows  past 
it  as...”  irrotational  flow  around  the  airfoU  plus  the  initial 


I 


■ 


i 


i 


I 


! 


i 


rolled  up  sheet  of  discontinuity.  This  avoids  the 
difficulty  of  infinite  velocities,  but  he  does  not  describe 
how  the  sheet  of  discontinuity  forms.  The  Russian 
authors  Kochin,  Kibei  and  Roze°  do  not  allude  to  the 
details  indicated  above,  but  only  say  that  "the  process  of 
shedding  of  vortices  with  prevdence  of  vortices  of  one 
sign  will  continue  until  the  circulation  f  reaches  the 
value  for  which  the  velocity  at  the  trailing  edge  of  the 
airfoil  will  be  finite".  These  differing  descriptions  are 
consistent  with  the  lack  of  direct  experimental  evidence 
about  the  origin  of  lift.  Up  to  now,  it  has  been 
impossible  to  measure  the  instantaneous  vorticity  over 
the  starting  vortex,  and  at  other  relevant  places.  To  our 
knowledge,  only  one  set  of  measurements  of  the 
circulation  in  the  fluid  downstream  of  an  impulsively 
started  airfoil  have  been  made.  Walker^'  used  a  dense 
concentration  of  particles  to  determine  the  growth  of  the 
total  circulation  in  the  fluid  aft  of  an  aiifoil  that  was 
started  from  rest  He  found  that  more  than  six  chord 
lengths  were  required  before  the  circulation  reached 
90%  of  the  steady  flow  value. 

2.  EXPERIMENTAL  TECHNIQUES 

2.1.  Review  of  the  mea.suremeni  technique 

The  technique  of  Laser  Induced  Photochemical 
Anemometry  (LIPA)  is  reviewed.  The  first  use  of 
excited  photochromic  chemicals  to  measure  fluid 
properties  was  made  by  Popvich  and  Hummel  They 
used  the  reaction  to  produce  time  lines,  from  which  they 
obtained  velocity  fields  for  some  quasi  one-dimensiomd 
flows.  The  extension  described  by  Falco  and  Chu^ 
which  is  reviewed  below  --  enables  time  resolved 
measurement  of  planar  velocity  vectors,  instantaneous 
Reynolds  stress  components,  the  in-plane  strain-rate 
components  and  the  component  of  vorticity 
perpendicular  to  the  plane  of  interest  in  fluctuating 
three-dimensional  velocity  fields. 

2.1.1.  The  photochromic  dve 

Use  is  made  of  the  phenomena  of 
’photochromism’  to  measure  both  velocity  and  vorticity 
in  the  starting  airfoil  problem.  Brown^  has  defined  it  as 
’a  reversible  change  of  a  single  chemical  species 
between  two  states  having  distinguishably  different 
absorption  spectra,  such  change  being  induced  in  at  least 
one  direction  by  the  aaion  of  electromagnetic 
radiation."  The  chemical  used  in  this  study  was  Kodak 
l,3-Trimethyl-8-nitrospiro(2-H-l-ben2opyran-2,2’- 
indoline}.  We  excited  it  at  351nm  using  pulses  of  20ns 
produced  by  a  Tachisto  Excimer  laser  running  on  XeF 
(approximately  lOOmj  output).  The  chemical,  which  was 
dissolved  in  kerosene,  reflected  blue  light  for  the  order 
of  tens  of  seconds.  Thus,  the  fluid  was  tagged  long 
enough  for  it  to  move  appreciably.  In  practice,  for  the 
velocities  of  our  experiments,  we  needed  to  follow  the 
lines  only  a  fraction  of  a  second  to  obtain  the 
information  required. 

2..L2,  Use  of  a  grid  to  tag  fluid  particles 


velocity  in  a  flow  using  the  visual  technique,  we  must  tag 
and  follow  specific  fluid  particles.  The  simplest  way  to 
do  this  is  to  have  two  Unies  intersect.  That  intersection 
defines  a  point  in  the  flow  which,  because  it  is  colored, 
we  can  follow.  To  measure  a  spatial  gradient,  we  must 
simultaneously  mark  two  points  in  the  flow.  To  enable 
us  to  measure  a  spati^  gradient  in  a  prescribed 
direction,  we  need  several  points  and  must  use 
interpolatioa  Thus,  the  creation  of  a  grid  of  intersecting 
marked  fluid  lines  enables  both  velocities  and  their 
gradients  to  be  measured. 

We  required  a  method  of  dividing  a  beam  of  laser 
light  into  ’n’  beams  where  ’n’  would  be  5-25.  In  doing  so, 
we  sought  to  lose  as  little  light  as  possible.  Thus,  we 
designed  a  divider  that  appears  as  an  oversized  bl^d 
reflection  grating  with  step  size  approximately  2mm  (see 
Figure  2).  The  incident  laser  energy  was  directed  at  a 
shallow  angle  to  the  plane  of  the  grating.  Because  the 
facets  are  so  much  larger  than  a  wavelength  of  light,  we 
expected  to  see  primarily  specular  reflection.  The 
separation  of  the  specularly  reflected  portions  of  the 
beam  were  determined  by  the  difference  of  the  angle  of 
incidence  and  the  blaze  angle,  and  the  grating  step  width 
and  spacing.  It  was  fabricated  by  individually  fixing 
aluminum  coated  mirrors  (with  the  coating  optimized  for 
shallow  incidence  angles)  to  a  machined  base.  Apart 
from  diffiaction  effects,  no  incident  laser  energy  would 
be  lost. 

As  expetted,  however,  diffiaaion  effects  were 
exhibited.  We  determined  that  as  long  as  the  grid  was 
formed  closer  than  approximately  1  m  from  the  divider, 
we  were  in  the  Fraunhofer  regime,  and  we  obtained  no 
noticeable  diffraction  effects. 


2.1.3.  Algorithms  to  obtain  fluid  dynamic 
quantities 

For  convenience,  we  pick  up  one  ’grid  box’  which, 
if  the  mesh  size  is  small  enough,  can  be  thought  of  as  a 
text  book  fluid  ’particle’  in  the  flow  at  time  t=0.  As  this 
particle  moves  with  the  flow,  it  may  undergo  translation, 
rotation  and  deformation.  Because  we  know  the  history 
of  this  specific  fluid  particle,  we  can  compute  the 
velocity  of  each  comer  by  taking  the  time  derivative  of 
displacement  ds  associated  with  each  point,  that  is, 
ds/Ai.  The  deduction  procedures  are  described  in  Falco 
and  Chu^.  In  a  three-dimensional  flow  field  a  fluid 
particle  may  have  motions  about  all  three  coordinate 
axes;  however,  we  will  obtain  the  projections  of  this 
motion  onto  the  plane  of  ihc  photos  which  is  parallel  to 
the  initial  plane  of  the  grid. 

To  obtain  velocity  gradients,  vorticity,  and  strain 
rates  from  the  velocities,  we  must  difference  across  the 
distance  of  a  grid  box.  Thus,  for  a  given  flow  the  grid 
must  be  small  enough  to  resolve  the  smallest  scales  of 
motion. 

An  alternative  method  can  be  used  if  we  are 
only  interested  in  the  vorticity.  Using  the  definition  of 
the  circulation  T . 

r  3  f  y.ds 


To  enable  us  to  measure  two  components  of  the 
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and  Gauss’  theorem  to  relate  the  surface  integral 
to  an  area  integral. 
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we  can  obtain  the  vorticity  component  normal  to 
the  film  plane,  from  a  grid  box  of  arbitrary  shape, 
with  circumference  C  and  area  A. 


The  velocity  component  V ‘ds  along  any  side  of  a 
grid  box  in  Eq.  (1)  is  estimated  by  forming  the  average 
of  the  velocity  components  obtained  at  two  of  the 
comers  of  the  box.  Dividing  F  by  the  area  of  the  box 
results  in  the  average  vorticity  at  the  centroid  of  this 
area.  This  technique  has  the  advantage  of  avoiding  a 
second  differencing  of  the  experimental  data. 

Thus,  by  following  the  distortions  of  a  single  grid 
box,  we  can  obtain  the  velocities,  circulation,  vorticity 
and  even  the  instantaneous  Reynolds  stresses  and  strain 
rates  if  desired,  over  that  small  area.  By  doing  this  for 
the  other  mesh  elements,  we  can  obtain  the  important 
fluid  dynamic  quantities  at  many  locations  over  a  two- 
dimensional  field  in  a  fluid  flow  simultaneously. 


2.2.  Experiments  performed 


Two  experiments  are  described.  The  first  may  be 
considered  as  a  calibration  of  the  UFA  technique.  Its 
objective  was  to  measure  the  boundary  layer  that 
develops  on  an  impulsively  started  moving  wall.  Since 
this  has  an  exact  solution  of  the  Navier-Stokes  equations, 
we  can  use  it  to  determine  the  accuracy  of  our  vorticity 
measurements.  The  second  experiment  concerned 
.leasurement  of  the  lift  formation  process  on  an 
impulsively  started  airfoil. 
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Figure  2.  The  optical  setup  used  to  produce  a  grid  of 
laser  lines. 


The  photographs  were  digitized  both  manually 
using  a  NAC  motion  analyzer,  and  automatically  using  a 
Westinghouse  ETV-2000  camera  which  was  controlled 
by  a  MegaVision  1024XM  image  processor.  The 
scanner  and  processor  both  are  capable  of  dealing  with 
images  as  large  as  2048  x  2048  pixels.  The  analysis 
algorithm  used  the  fact  that  relatively  straight  lines  were 
present  in  the  flow,  and  the  intersections  of  those  lines 
were  the  points  to  be  followed.  By  using  spatially  biased 
filtering  kernels,  the  centerline  of  each  set  of  lines  was 
extracted  from  the  image,  and  the  logical  interse«ion 
".AND"  was  obtained. 

The  intersection  information  was  processed  into 
fluid  kinematic  quantities  on  an  IRIS  3120.  The  IRIS 
provided  not  only  a  fast  CPU  but  also  a  graphics 
environment  in  which  the  individual  steps  of  the 
algorithm  could  be  checked  for  consistency. 


2.1.4.  Error  analysis 


Consider  a  single  grid  box.  We  will  assume  that 
the  grid  mesh  is  spaced  e  and  that  a  worse  case  line 
movement  between  images  is  1/2  the  mesh  width.  After 
time  At  standard  analysis  (see  Taylor  19821  shows  that 
the  error  in  estimating  the  vorticity  is 

So,^  =  (l/M)(5t/e)  (3) 


to  first  order  in  e . 

As  an  example  of  the  designed  accuracy  of  our 
experiments,  if  we  produce  lOO^m  lines  on  a  grid  mesh 
of  1mm,  and  our  error  in  reading  the  center  of  the  lines 
is  10%.  we  obtain  6«/e  =  .01.  For  At  =  .Olsec,  the 
error  in  the  measurement  of  vorticity  is  1/sec.  This  is 
better  than  has  been  possible  with  hot-wire  probes,  and 
is  comparable  with  the  most  sophisticated  single  point 
laser  doppier  measurements  (see  Lang^*^). 


For  both  experiments  the  same  laser,  chemicals 
and  similar  optical  arrangements  were  used.  A  Tachisto 
excimer  laser  outputting  lOOmj  and  351nm  was  used  to 
excite  the  chemi*^.  Figure  2  shows  a  sketch  of  the 
optical  semp  used  to  produce  a  grid  of  laser  lines. 
Figure  3  shows  an  example  of  the  grid  produced  after  it 
has  been  distorted  by  the  fluid  motion.  For  the  Stokes’ 
layer  problem  the  grid  extended  across  the  entire  layer. 
For  the  airfoil  problem  Figure  4  shows  the  size  of  the 
measurement  volume  relative  to  the  starting  vortex 
diameter  as  indicated  by  the  roll  up  of  photochromic  dye 
generated  on  the  airfoil  for  the  case  when  the  airfoil  had 
traveled  1.3  chords.  This  is  the  largest  it  grew  to  during 
our  measurements.  Thus  we  anticipate  capturing 
essentially  all  of  the  vorticity  in  the  starting  vortex. 
Figure  5  shows  the  setup  used  to  create  the  impulsively 
started  flows.  The  mylar  belt  was  chain  driven  by  a  1/4 
HP  Dayton  DC  motor  regulated  by  a  Dayton  SRC 
controller.  In  the  starting  airfoil  experiment  this 
combination  accelerated  the  belt  to  a  final  velocity  of 
61mm/s  in  .06  sec  (.057  of  a  chord  length).  The  airfoil, 
which  approximated  a  NACA  0012,  had  a  chord  length 
of  63.5mm.  The  chord  Reynolds  number  was  1550,  and 
it  was  fixed  at  a  4.7  degree  angle  of  attack. 
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A  35mm  high  speed  Photosonics  35ML  camera 
with  a  quartz  regulated  motor  was  used  to  record  data. 
A  lens  arrangement  providing  approximately  1:1 
magnification  was  used.  For  the  Stokes’  flow 
experiments,  we  framed  at  150  frames/second,  while  for 
the  airfoil  experiment  50  frames/second  were  used.  For 
the  airfoil  problem  the  following  timing  sequence  was 
used.  The  camera  was  started  approximately  1  second 
before  the  belt  to  bring  it  up  to  speed,  then  the  belt 
motor  was  started  along  with  the  delay  timer.  The  delay 
timer  triggered  the  laser  at  preselected  delay  times 
corresponding  to  .3,  .8  and  1.3  chord  lengths  of  travel  of 
the  airfoil.  Only  the  first  two  frames  after  the  laser  fired 
were  used,  in  each  case,  for  the  calculation  of  the 
vonicity. 
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We  can  test  the  accuracy  of  the  technique  in  a 
flow  for  which  an  exact  solution  of  the  Navier-Stokes 
equations  is  obtainable.  The  flow  is  that  of  an 
impulsively  started  infinite  plate  in  an  initially  still 


Figure  3.  A  photo  of  the  grid  after  it  has  been  distorted 
by  fluid  motion.  The  wing  is  moving  left  to  right. 


Figure  4.  The  size  of  the  measurement  region  with 
respect  to  that  of  the  starting  vortex  at  .3,  'S  and  1.3 
chords  behind  the  airfoil.  Tlie  wing  is  moving  left  to 
right. 


Figure  3.  Setup  used  to  create  and  measure  the  Stokes' 
ilow.  This  same  device  was  modified  to  support  an 
airfoil  (see  inset). 


infinite  ambient.  The  problem  involves  the  two- 
dimensional  diffusion  of  vorticity  generated  at  the  plate 
into  the  ambient.  The  thickness  of  the  boundary  layer, 
S,isa  function  of  the  kinematic  viscosity,  i/,  and  the  time 
after  the  start  of  the  plate’s  motion.  We  have 
approximated  this  motion  in  a  tank  using  a  belt  which 
can  be  rapidly  staned  (see  Figure  5).  The  width  of  the 
belt  is  17.8  cm,  and  the  distance  between  the  rollers  is 
152.4  cm,  with  the  test  position  90  cm  downsueam  of  the 
lead  roller.  All  experimental  data  was  obtained  before 
the  leading  edge  effects  reached  the  test  position. 
Furthermore,  the  aspect  ratio  is  sufficient  to  prevent  the 
disturbances  generated  in  the  upstream  comers  from 
reaching  the  center  of  the  belt  at  the  test  section  in  the 
test  time  interval.  Therefore,  the  wall  layer  flow 
measured  was  two-dimensional.  The  speed  of  the  belt 
was  adjustable  from  2.5  to  22.9  cm/sec.  It  reaches  a 
constant  speed  within  .06  second.  This  short 
acceleration  reasonably  approximates  an  impulsive  start. 
A  grid  mesh  size  of  .155  was  used,  with  the  laser  lines 
intersecting  as  in  Figure  5. 

Results  of  the  velocity  gradient  obtained  from 
double  differentiation  in  the  Stokes’  layer  are  shown  in 
Figure  6.  where  they  are  compared  with  the  exact 
solution.  Both  data  reduction  by  hand  and  using  the 
image  processing  equipment  are  shown.  All  quantities 
are  non-dimensionalized  by  the  similarity  variables  of 
the  exact  solution.  Because  the  flow  is  two-dimensional, 
incompressible,  and  doesn’t  change  with  the  streamwise 
coordinate,  x,  the  continuity  equation  reduces  to  dv/dy 
=  0,  and  a  check  on  this  is  also  shown  in  Figure  6  (v  is 
the  velocity  component  normal  to  the  wail).  We  see  that 
both  the  error  in  the  estimates  of  the  continuity 
equation,  and  in  the  velocity  gradient  distribution  are  of 
similar  magnitude. 
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Figure  b.  Comparison  of  experimental  measurements  of 
the  velocity  gradient  in  a  Stokes’  layer  compared  with 
the  theoretical  solution.  (  +  )  the  algorithm,  (o)  manually. 
(-)  the  exact  solution.  A  check  of  the  continuity  equation 
dv/dy  =  0  is  also  shown,  (x)  the  algorithm,  (' j  manually. 


We  can  make  use  of  all  of  the  samples  in  the 
continuity  equation  to  get  a  statistical  estimate  of  the 
eiTor  in  measurement  of  velocity  gradients.  The  rms 
error  equals  .9/sec,  consistent  with  the  design  error 
analysis.  The  deviation  of  the  experimental  profile  in 
the  near  wall  region  is  a  consequence  of  the  grid  mesh 
being  too  large  to  adequately  resolve  the  large  value  of 
the  gradient  du/dy.  On  the  whole  the  results  indicate 
that  we  should  be  able  to  measure  vorticity  over  an  area 
more  accurately  than  instruments  using  multiple  hot¬ 
wires  or  multiple  LDA  beams  can  at  a  point. 
Furthermore,  the  resolution  of  each  vorticity 
measurement  can  potentially  be  higher,  since  we  are 
essentially  constrained  only  by  the  minimum  thickness  of 
an  ultraviolet  laser  beam. 

3.2.  Impulsively  started  airfoil 

Using  the  photochromic  dye,  we  could  observe 
the  direction  of  motion  of  fluid  particles,  and  measure 
the  time  required  for  a  point  in  the  fluid  to  move  along 
the  upper  surface  to  the  trailing  edge.  We  could  observe 
motion  in  increments  of  Ut/c  =  .02.  The  motion  of 
marked  fiuid  was  always  observed  to  be  downstream 
with  respect  to  the  airfoil.  By  positioning  the  airfoil 
under  the  laser  generated  grid,  and  firing  it  at  a 
repetition  rate  high  enough,  we  could  successively  mark 
selected  "fluid  particles"  on  the  airfoil’s  surface  and 
follow  their  movements.  Doing  this,  we  obtained  a 
definite  indication  that  the  separating  sueamline  was 
fixed  to  the  trailing  edge  by  Ut/c  =  .06.  We  could 
resolve  movement  of  .006  chord.  During  this  time  a 
vortex  formed  on  the  LOWER  surface  of  the  airfoil,  very 
close  to  the  trailing  edge,  but  beneath  the  fixed 
separation  point.  It  should  be  noted  that  although  the 
trailing  edge  was  sharp  when  viewed  with  respect  to  the 
airfoil  as  a  whole,  it  had  a  radius  of  curvature  of 
approximately  ,42mm. 

Having  established  that  the  process  of  movement 
of  the  separation  point  was  completed  in  Ut/c  =  .06,  we 
examined  the  starting  vortex.  Using  the  grid  technique 
described  above,  we  measured  the  vorticity  and 
circulation  in  the  starting  vortex  after  it  had  convected 
.3,  .8,  and  1.3  chord  lengths  downstream,  r/r*,  was 
.392,  .378  and  J87.  Thus,  the  magnimde  of  the 
circulation  has  not  changed  within  our  experimental 
error.  The  steady  state  circulation  fpt  our  NACA  0012 
airfoil  was  estimated  to  be  1024mm“'  sec.  The  estimate 
at  the  .3  chord  included  all  of  the  circulation  aft  of  the 
airfoil.  We  see  that  the  estimates  indicate  the  staning 
vortex  has  developed  to  full  strength  by  .3  chord.  Thus, 
the  circulation  in  the  starting  vortex  is  approximately  1/3 
of  the  total  circulation  around  the  airfoil  in  steady 
motion. 

The  fact  that  r  of  the  startine  vortex  had  not 
reached  the  theoretical  value  of  xcU^sina  by  this 
position  downstream  was  expected  from  residts  of 
Walker^  ^  and  theoretical  results  of  Wagner^  .  and 
numerical  calculations  of  Basu  and  Hancock^ and  Chow 
and  Huang  ,  among  others.  Explanations  have  alluded 
to  the  fact  that  when  the  starting  vortex  was  at  these 
chord  positions,  the  separation  streamline  had  not  yet 
reached  the  trailing  edge,  so  that  additional  vonicity  was 
being  shed  ~  which  might  rearrange  itself  into  additional 


vortices.  However,  we  find  that  the  separating 
streamline  has  moved  to  the  trailing  edge  and  remains 
fixed  there  by  .06  chords,  so  the  classical  mephanitin  for 
the  supply  of  vorticity  is  no  longer  operative. 

Figure  7  compares  our  results  with  those  of 
experiments  of  Walker^'  and  theoretical  estimates  of 
Wagner^°.  It  is  interesting  to  note  that  the  value  of  the 
circulation  we  measured  at  .3  chord  is  (within 
experimental  error)  the  same  as  that  measured  by 
Walker  and  calculated  by  Wagner.  They  are  at  all  rimes 
considering  the  total  cir^ation  of  the  fluid  downstream 
of  the  airfoils.  At  3  chord  we  are  also  measuring  the 
total  circulation  of  our  airfoil.  This  correspondence  is 
expected,  and  adds  weight  to  the  accuracy  of  our 
measurements.  The  discrepancy  found  at  .8  and  13 
chords  thus  highlights  the  fact  that  processes  other  than 
the  classical  physical  piemre  are  taking  place.  It  is 
interesting  to  point  out  that  for  impulsively  started 
airfoils,  only  a  single  vortex  is  expected  if  sii^arity  is 
assumed  (Pullin^^).  Thus,  it  appears  that  the  process  of 
roll  up  of  the  vorticity  shed  from  the  back  of  the  airfoil 
does  not  conform  to  this  hypothesis,  since  the  starting 
vortex  contains  only  a  third  of  the  vorticity  responsible 
for  the  lift  of  the  airfoil. 

Further  experiments  were  performed  to 
determine  whether  additional  vorticity  was  being  shed 
off  the  airfoil  after  the  circulation  in  the  starting  vortex 
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Figure  7.  Tire  increase  in  circulation  with  distance  from 
the  airfoil,  compared  with  experimental  and  theoretical 
data. 
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bad  reached  its  constant  value.  A  single  photochromic 
line  was  created  .07  chord  lengths  behind  the  airfoil  after 
it  had  traveled  .5,  .75,  1  and  2  chord  lengths.  We 
examined  the  movement  of  this  line,  both  above  and 
below  the  airfoil  to  determine  whether  a  shear  existed 
across  the  wake.  Indeed  a  shear  existed  in  ail  the 
measurements.  Figure  8  shows  the  distortion  of  a  single 
photochromic  line  created  approximately  .06  chord 
lengths  behind  the  airfoil,  and  photographed  after 
approximately  .22  chord  lengths  of  travel,  for  the  case 
where  the  airfoil  has  traveled  1/2  chord,  i.e.  after  the 
starting  vortex  has  attained  its  final  value  (see  Figure  7). 
The  sign  of  this  vorticity  was  the  same  as  that  in  the 
starting  vortex.  This  magnitude  decreased  as  the  airfoil 
had  traveled  further,  but  measurable  shear  existed  in  the 
experiment  even  after  two  chords  of  travel.  Thus,  there 
was  continued  generation  of  vorticity  long  after  the 
separation  streamline  was  fixed  to  the  trailing  edge,  and 
the  Kutta  condition  was  therefore  not  satisfied  after  at 
least  2  chords  of  travel  in  our  measurements. 


Figure  8.  The  shear  that  exists  across  the  wake  when  the 
starting  vortex  is  0.5  chord  downstream.  The  wing  is 
moving  right  to  left. 


4.  DISCUSSION 


The  demonstrated  accuracy  of  the  technique  is 
sufficiently  adequate  to  enable  us  to  reexamine  many 
unsteady  fluid  dynamic  problems.  Since,  in  principle, 
the  calibration  depends  only  on  a  time  and  a  length 
scale,  there  is  further  room  for  improvement.  In 
practice,  determining  the  intersection  points  of  the  grid 
lines  is  limited  by  the  grain  size  of  the  film  used  (or  the 
resolution  of  a  photo  optical  detector),  the  distortion  of 
the  developing  process  and  the  accuracy  of  either  hand 
or  automated  measurements  in  defining  the  lines. 
Improvements  in  the  grain  limitation  can  be  made  by 
using  a  larger  format  camera  or  larger  aperture 
lens/finer  grain  film  combination. 

For  these  preliminary  low  Reynolds  number 
experiments,  the  average  magnitude  of  the  vorticity  in 
each  of  our  grid  boxes  was  3  sec"  .  Since  the  accuracy  of 


our  technique  is  2  .9  sec'^  we  could  only  expect  66% 
accuracy  per  individual  estimate.  However,  estimates  of 
the  circulation  for  each  case  involved  integrating  over  an 
average  of  23  circulation  cells,  thus  significantly  reducing 
the  error  in  estimating  r.  Furthermore,  in  a  direct  test 
of  the  reproducibility  of  our  estimates,  at  .8  chords  two 
independent  runs  were  analyzed.  These  results 
indicated  that  the  circulation  varied  by  4%.  Another 
possible  source  of  error  was  the  possibility  of  not 
including  all  of  the  vorticity  in  the  staning  vortex  for  the 
case  at  1.3  chords.  However,  as  Figure  4  indicates,  the 
dye  suggests  that  most  of  the  vorticity  has  been 
measured. 

4.1,  The  generation  of  lift 

If  the  circulation  is  no  longer  determined  by  the 
fact  that  the  separating  streamline  has  moved  to  the 
trailing  edge,  then  what  determines  the  magnitude  of  the 
circulation  responsible  for  steady  lift?  Where  is  the 
additional  circulation  coming  from?  In  the  theoretical 
and  numerical  accounts,  vorticity  is  continuously  coming 
off  in  a  sheet  that  rolls  up  into  the  starting  vortex  at 
oome  position  downstream  of  the  airfoil;  the  vorticity 
generation  comes  to  an  end  when  the  Kutta  condition  is 
established.  Graham'  critiques  analysis  which  suggest 
that  a  sheet  is  shed  approximately  parallel  to  the 
bisector  of  the  trailing  edge  and  then  rolls  up.  He  notes 
that  all  evidence  suggests  that  the  surface  of 
discontinuity  rolls  up  immediately  as  it  forms  and  comes 
off  the  airfoil.  Using  this  hypothesis  he  arrives  at  a 
picture  of  the  initia'  lift  being  very  high,  decreasing  and 
then  increasing  as  the  Wagner^®  function  after  Ut/c  = 
.1.  All  of  this  work  suggests  that  as  the  circulation  is 
increasing  the  rear  stagnation  point  is  moving  towards 
the  trailing  edge.  However,  it  takes  several  chord 
lengths  to  have  all  of  the  vorticity  shed. 

The  fact  that  we  measured  only  2/3  of  the 
circulation  experted  in  the  starting  vortex  after  1.3 
chords,  and  that  there  does  not  appear  to  be  any  change 
in  the  strength  over  the  previous  chord  of  motion, 
suggests  that  additional  vorticity  is  continuing  to  be  shed 
after  the  vortex  rolls  up.  A  contradiction  arises.  The 
classical  theory  suggests  that  the  upper  stagnation  point 
will  have  moved  to  the  trailing  edge  when  the  starting 
vortex  has  been  shed,  and  thus  there  is  no  longer  a 
mechanism  for  the  generation  of  additional  circulation. 
The  unique  irrotational  solution  for  the  circulation  is 
established.  However,  in  the  irrotational  solution  the 
separating  streamline  is  a  stagnation  streamline;  in  the 
real  fluid  it  is  not. 

In  our  experiments  the  separation  point  has 
moved  to  the  trailing  edge  after  .06  of  a  chord  of  travel. 
Thus,  it  appears  that  the  movement  of  the  separation 
point  to  the  trailing  edge  does  not  in  itself  establish  the 
Kutta  condition.  It  must  be  cautioned  that  we  are  not 
able  to  resolve  the  velocity  field  with  sufficient 
resolution  to  determine  the  stagnation  point,  but  our 
determination  is  being  made  by  measurirc  the 
movement  of  lines  of  marked  fluid  that  are  created  at 
large  angles  to  the  surface,  and  selected  marked  'fluid 
particles’  (intersections  of  the  marked  fluid  lines).  We 
see  that  initially  the  lines  and  marked  points  move 
downstream  with  respect  to  the  airfoil,  and  that  after  .06 
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chord  of  airfoil  travel  they  are  fixed  to  the  trailing  edge. 
Thus,  it  is  not  at  all  certain  that  the  initial  position  of  the 
stagnation  point  is  on  the  upper  wing  surface  (we  have 
no  evidence  for  or  against  this  possibility).  After 
movement  to  the  trailing  edge  (if  this  happens  at  all), 
one  still  has  the  presence  of  a  strong  shear  layer. 
Accordingly,  we  feel  that  the  shear  layer  still  exists  for 
more  than  6  chords,  Le.  the  circulation  is  changing  over 
airfoil  travel  of  more  than  50  times  as  long  as  indicated 
by  the  classical  theory. 

Finally,  we  note  that  this  is  a  very  low  Reynolds 
number  experiment  It  could  be  argued  that  the  lift  on 
this  airfoil  will  not  be  the  same  as  that  for  a  NACA  0012 
at  higher  Reynolds  numbers.  However,  it  is  interesting 
to  note  that  the  value  for  the  total  circulation  which  both 
Walker  and  we  have  measured  after  .3  chords  is 
remarkably  close.  It  is  certainly  possible  to  perform 
experiments  of  this  type  at  considerably  higher  Reynolds 
numbers.  The  limitation  is  not  the  LIPA  technique,  but 
instead  the  mechanical  mechanisms  for  rapid 
acceleration  of  the  airfoil. 


5.  CONCLUSIONS 

The  generation  of  lift  on  an  airfoil  is  re-examined. 
Using  a  new  quantitative  optical  measuring  technique, 
we  determined  that  the  separation  streamline  moves  to 
the  trailing  edge  of  an  airfoil  long  before  the  Kutta 
condition  (unsheared  flow  from  the  trailing  edge)  is 
established,  i.e.  long  before  the  change  in  circulation 
necessary  to  reach  the  steady  flow  lift  has  occurred. 
Preliminaiy  results  show  that  the  vorticity  in  the  starting 
vortex  which  has  convected  approximately  1.24  chords 
after  the  airfoil  separation  streamline  has  moved  to  the 
trailing  edge,  is  much  less  than  that  predicted. 

This  suggests  that  the  classical  hypothesis  relating 
the  formation  of  lift  and  the  establishment  of  the  Kutta 
condition  needs  to  be  re-examined. 
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Appendix  D  contains  the  orifainal  List  of  equipment  requested 


under  AFOSR-86-0242.  It  also  contains  the  list  of  equipment  which  was 
ultimately  purchased,  as  well  as  a  copy  of  the  letter  explaining  our 
reasons  for  requesting  equipment  substitutions. 


71-2309 

MSU  EQPT. 

COST -SHARING 

MSU  OTHER 
COST- SHARING 

CAM  Corp.  (P.O.  36154) 

Seiko  color  copier  with 
video  multiplexer ; maintenance 

Helix  Ltd.  (P.O.  37409) 

Hasselblad  135MM  Macro  lens  for 
Westinghouse  Camera  System 

$  1,039.00 

$  9,741.00 

$  1,593.00 

1,053.00 

Helix  Ltd.  (P.O.  52681) 

Hasselblad  lens  mounting  ring; 
gelatin  filter  holder;  lens 
mount  adapter;  auto  bellows 
extension;  shipping 

Megavision  (P.O.  36153) 

Megavision  Image  Processing 
System,  1024XM;Dual  output 
channel ; 2 -4Mbyte  image 
memory ; Pixel  Data  Flow 

Processor; 2  1  Mbyte  boards, 

Color  monitor ; maintenance 

53,489.03 

3,008.41 

2,397.30 

Optimem  (P.O.  36155) 

Optimem  1000  Disc  Driver; 

4  singlesided  media; 
ext .  warranty 

13,600.00 

1,646.00 

Silicon  Graphics  (P.O.  36156) 

Iris  3120  workstation; 
tape  drive;  accelerator; 
processor ; compiler ; 
extra  display  memory; 
sof twar e ; ext .  warranty 

44,525.02 

5,590.01 

Sun  Microsystems  (P.O.  36152) 

Sun  3/260HM-P2  Workstation; 

481A  Multiplexer ; License , 

Sof tware ; Disc  Driver ;DNA/ 

2-01  Sof tware ; maintenance 

44,919.22 

2,900.00 

5,437.50 

Vanguard  (P.O.  36417) 

M-35-C  Projection  head  -  8,863.29 

Westinghouse  (P.O.  36416) 

ETV-2000  Camera  System; 

maintenance  19,800.00  1,000.00 


Travel  .  .  3,646.62 

Freight 


Other 


1,044.  50 
2,704.93 


Overhead  (42X  MTDC  thru  6/30/88; 
44X  7/1  thru  7/31/88) 


9.444.29 


TOTALS 


$  177.372.27 


$  28,211.70 


$  31,858.15 
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ORIGINAL  EQUIPMENT  LIST 


EQUIPMENT 


OIPIX  imaae  processinq  system 

Arles  III  image  processor  31,500 
3  megabytes  image  memory  15,600 
pixel  processor  9,950 
vector  graphics  package  10,000 
extended  software  package  15,000 
Hitachi  19"  color  monitor  5,200 
16-32  bit  image  converter  3,000 
fortran  77  2,500 
vgtlOO  graphics  monitor/printer  3,500 
installation  3,000 
Training  2,000 
support  i  year  at  1%/mo  11.550 


Subtotal  112.800 


EIKONIC  DIGITIZER 

EC78/99  Digital  camera  16.500 
L-01  Condenser  light  source  5,500 
C-10  Workstation  4,800 
E-12  interface  to  Microvax  II  2,500 
Driver  1,500 
EO-2  Data  Normalizer  1,000 
I-Ol  Color  filter  wheel  2,500 
previewer  1.800 
Disc  application  software  1,000 
105  mm  Nikon  macro  lens  900 


38.000 

Less  6%  educational  discount  -2,280 


35.720 

Factory  installation  1.800 

Support  1  year  at  1%/mo  4,344 


Subtotal  41.864 

DEC  MICRO  VAX  II  HOST 

Micro  Vax  II  32.686 


Includes : 

MVII  SBB:1MB/FP  BA123 
RD/RX  DISC  CONTROLLER 
RX50  W/ CABLES 
RD53  W/ CABLES 
MICROVMS  SOFTWARE 
NETWORKING 
MANUELS 
FORTRAN  77 


DEC  RA81  456  mb  disc /controller  27,000 

DEC  9-track  tape  drive  9,990 

69,676 

Less  14%  educational  discount  -9,755 

59,921 

8mb  Chryslln  memory  3,380 

installation  2,000 

support  1  year  at  1%/mo  7,190 

Subtotal  72,491 


Total  Equipment  202.071 

Total  Direct  Costs  227,155 

Indirect  Costs  (41%  of  Direct  Costs 

less  equipment)  10,284 

Total  Cost  237,439 


MSU  Cost  Sharing  (25.3%) 
Total  Cost  to  DOD 


-60,000 

$177,439 


MICHIGAN  STATE  UNIVERSITY 


OFFICE  OF  ASSOCIATE  DEAN  FOR  EAST  lANSING  •  MICFIJGAN  •  48M4-I226 

GRADUATE  STUDIES  AND  RESEARCH 
COLLEGE  OF  ENGINEERING 


February  2,  1987 


Dr.  Jim  McMichael 
AFOSR/NA 
Bolling  AFB 

Washington,  D.  C.  20332-6448 
Dear  Jim: 


Request  for  Equipment  Substitutions  in  POD  Equipment  Grant 

rtFOSR-86-0242  DEF 


I  have  made  a  thorough  re-evaluation  of  the  image  processing  equipment  in 
the  grant,  and  have  determined  that  a  better  system  can  be  put  together  for 
the  same  costs.  This  system  will  enable  us  to  more  accurately  measure  the 
vorticity  and  strain  rates,  and  to  better  understand  the  data.  It  will  also 
enable  us  to  interface  with  the  Direct  Numerical  Simulations  that  will  be 
coming  from  NASA  Ames  and  other  places. 

The  major  new  components  of  the  system  are: 

a.  Westinghouse  2"  video  camera  that  gives  2048  x  2500  pixels; 

b.  Megavision  1024M  processor  that  will  process  2048  x  2048  images  and 
can  be  upgraded  to  4096  x  4096  images  which  will  be  necessitated  by 
later  higher  speed  experiments; 

c.  Silicon  Graphics  IRIS  processor,  which  will  allow  us  to  view  the 
data  from  any  angle  and  to  reconstruct  the  three-dimensional  fields 
from  the  two-dimensional  slices  of  vorticity  data  we  are  obtaining, 
and  compare  them  with  numerical  computations. 

d.  Sun  3/260  with  optical  disk,  to  be  the  host  computer  for  the 
network  of  digitizer,  image  processor  and  graphics  processor,  mass 
storage,  and  communications. 

e.  A  hard  copy  device  that  will  enable  permanent  records  at  the 
resolution  of  the  screen  images. 

A  detailed  list  of  the  new  equipment  is  attached.  In  each  case  I  have 
bargained  for  some  special  deal  with  the  vendors.  The  net  result  is  that  I 
have  achieved  this  significant  enhancement  of  system  capabilities  at  no 
additional  cost. 
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Permission  is  requested  to  make  these  substitutions. 
Sincerely, 


Robert  Falco 
Professor  and  Director 
Turbulence  Structure  Laboratory 

Approved : 


John  R.  Lloyd,  Chairperson  Howard  G.  Glider,  Director 

Mechanical  Engineering  Dept.  Contract  &  Grant  Administration 


William  C.  Taylor,  Acting  Assoc.  Dean 
Graduate  Studies  and  Research 
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PROPOSED  REVISED  EQUIPMENT  COST  BREAKDOWN 
FOR  AFOSR-S6-0242,  R.  E.  FALCO 
1/26/87 


TOTAL  COST 


GRANT  MSU  EQ  MSU  S&S 
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